This notebook was created in Jupyter lab. The required Python packages are stored in the environment file environment_gastrushape.yml. To import this environment, see the instructions in step 0.1 and 0.2, but replace the name of the yml-file with environment_gastrushape.yml and the environment name with gastrushape.
# The basic libraries
import numpy as np
import pandas as pd
import os
import math
import scipy
# To plot data
import matplotlib
from matplotlib import pyplot as plt
import matplotlib.image as mpimg
matplotlib.use('Agg')
import sys
# To read in images
import imageio
# To process images
import cv2
# To read ROIs generated in ImageJ
from read_roi import read_roi_file
from read_roi import read_roi_zip
# To fit ellipse
from skimage.measure import EllipseModel
from matplotlib.patches import Ellipse
We obtained the LOCO-EFA values from gastruloid shapes that were extracted from Nikon ND2 files. Our ND2 files often contain images from multiple gastruloids in one field of view, of which we recorded multiple fluorescent channels, multiple positions in Z (z-stack), and multiple positions in XY. We developed a segmentation tool that deals with this format and extracts outlines from individual gastruloid shapes.
To run this tool, follow below steps. The tool may also be used for image files with a different extension, such as .jpg, .png, .tif or .tiff. If the image file contains a bright-field image, segmentation will likely not be performed correctly. The signal in the image needs to be limited to the gastruloid and should not be visible in the background. Therefore, the image file should only contain fluorescent channels, ideally one that is homogeneously showing the gastruloids shape, for example by using DAPI to visualize all nuclei.
NB: Make sure that your images contain unique gastruloids, or that you remove 'doubles' from the analysis later on. In case you imaged multiple gastruloids and took several images with overlapping areas, there's a chance that two different images include the same gastruloid. Make sure you remove 'doubles' from the LOCO-EFA analysis to not have them in the final plot. We carefully checked for 'doubles' and did not include them in the plot.
All the required files are stored in the Interface folder. Navigate in your terminal to this folder using the cd tool. In this folder you will find the following files:
How to run the segmentation tool gui-segmentation.py?
You first have to create the environment from the environment.yml file by running
conda env create -f environment.yml
The name of the environment is ciha. You can check if it's in the conda environments list by running
conda env list
You then activate the environment by running
source activate ciha
The next time you want to use the segmentation tool, you only need to navigate to the interface folder and run step 2 and 3.
Explanation of buttons (in order of use)
File Open: opens a file dialog window to let you select a folder. If the folder contains multiple image files, all will be recognized when they have the extension .nd2, .png, .jpg, .tif, or .tiff.
Run segmentation: runs the segmentation on the first image in the folder with the standard parameters (threshold value = none, minimal distance = 70, minimal size = 100, Gaussian smoothing = 5, opening for maxima = 5 [pixels?]). These parameters can be adjusted if desired.
For each image position (XY), save the outlines:
File Save Outlines: opens a file dialog window. Select the folder in which the outlines should be stored and provide a name for the file including the extension (e.g., .txt or .csv).
Save Outlines: saves the outlines that were found after segmentation in the selected folder with the set name.
Next xy: If an image file contains multiple XY positions, you may click Next xy to go to the next position and repeat the segmentation procedure, the selection of the folder and naming the file, and saving the outlines.
Next File: if the folder contains multiple images, you select the next image by clicking this button. The segmentation tool keeps track of the image number and the position number in the boxes Current Image and Current xy, respectively.
Run Fourier Analysis: runs a Fourier analysis on the outlines. The calculated Fourier components can be stored using the buttons File Save Fourier to select the path and file name, and Save Fourier Components to store the components in the selected path.
The segmentation tool stores all the outlines of the detected gastruloids in a single .csv or .txt file. To generate from outlines LOCO-EFA values that describe the shape, each outline should be stored as a separate .csv file. To do so, you may go over the following part of this jupyter notebook that stores the outlines separately.
# Import the required libraries
from numpy import array, savetxt
# Select the txt file (a file dialog will open)
try:
from tkinter import Tk
from tkFileDialog import askopenfilenames
except:
from tkinter import Tk
from tkinter import filedialog
Tk().withdraw() # we don't want a full GUI, so keep the root window from appearing
filenames = filedialog.askopenfilenames() # show an "Open" dialog box and return the path to the selected file
filename = ''.join(filenames)
print(filename)
# Read in file, and get the number of outlines listed in the file
reader = open(filename,'r')
try:
#print(reader.read())
value_lines = list(reader)
finally:
reader.close()
size_val = len(value_lines)
# Get the number of outlines listed in the file and check, if not correct, adjust it
vals = value_lines[size_val-1].strip("\n")
vals2 = vals.split("\t")
n_outlines = int(vals2[0][-1])+1
print('Number of outlines in this file are:', n_outlines)
#n_outlines = 3
#print(n_outlines)
# Create a new folder to store the separate outlines in
# Get the file path (file_path) and file name (report_name) without the extension
filename_split = filename.split("/")
filename_split_length = len(filename_split)
file_path = (filename[0:(len(filename)-len(filename_split[filename_split_length-1]))])
print(file_path)
report_nametxt = filename_split[filename_split_length-1]
report_name = report_nametxt[0:(len(report_nametxt)-4)]
print(report_name)
# Create a new folder to store the header and the separate outline files in
newpath = '%sFolder_%s' % (file_path,report_name)
print(newpath)
if not os.path.exists(newpath):
os.makedirs(newpath)
In case you have many imaging files, you will also have many gastruloid outlines. To save outlines with a file number starting from a certain number (to eventually have outline_1 to outline_100 instead of 10 outline_1 files), change count_start to the desired number. Keep a record on this starting number for each new XY, and subsequently each new imaging file. Note that the code in the cell below doesn't work yet for txt or csv files that contain more than 10 outlines. To handle those, first split the files in two (0-9 and 10-19) or more (20-29, etc.).
>> currently the code only reads out the last value of the csv or txt file (numbers 0-9) to determine the number of outlines that are listed in this file.
# To save outlines with a file number, starting from a certain number, change count_start to the desired number.
count_start = 10;
# For each outline (j), list the x and y values in separate arrays and write the arrays to a csv file
for j in range(0,n_outlines):
x_array = np.array([])
y_array = np.array([])
for i in range(1,size_val):
# print('i= %i and j= %j',(i,j))
value_lines_stripped = value_lines[i].strip("\n")
value_lines_split = value_lines_stripped.split("\t")
if j == (n_outlines-1) and i == (size_val-1): # j is the last outline and i is the last written line
# >> append array, save array, save header
xval = float(value_lines_split[1])
xvalr = round(xval)
yval = float(value_lines_split[2])
yvalr = round(yval)
x_array = np.append(x_array,xvalr)
y_array = np.append(y_array,yvalr)
header_outline = value_lines[0].strip("\n")
np.savetxt('%s/Outline_%i.csv' % (newpath, j+count_start), np.array([x_array, y_array]).T, fmt='%i', delimiter=',')#, header = header_cont)
np.savetxt('%s/Header_%s.csv' % (newpath,report_name), [0], header = header_outline)
elif int(value_lines_split[0][-1]) == j: # if j is not yet 9, and line i belongs to j >> append array
xval = float(value_lines_split[1])
xvalr = round(xval)
yval = float(value_lines_split[2])
yvalr = round(yval)
x_array = np.append(x_array,xvalr)
y_array = np.append(y_array,yvalr)
elif int(value_lines_split[0][-1]) > j: # if j is not yet 9, and line i is larger than j >> save text
np.savetxt('%s/Outline_%i.csv' % (newpath, j+count_start), np.array([x_array, y_array]).T, fmt='%i', delimiter=',')#, header = header_cont
i = size_val
It can be helpful to visualize the gastruloid shape from the outline file. The following code allows for that.
outline_data = pd.read_csv('Outlines/outlines138', header=None)
n_pixels_outline = round(outline_data.size / 2)
outline_matrix = np.zeros((1024,1024))
for i in range(n_pixels_outline):
read_x = round(outline_data.loc[i,0])
read_y = round(outline_data.loc[i,1])
outline_matrix[read_x,read_y] = 1
# Fill the outline image
filled_binary = scipy.ndimage.binary_fill_holes(outline_matrix)
# Plot shape
%matplotlib inline
plt.imshow(filled_binary, cmap = 'gray')
<matplotlib.image.AxesImage at 0x7f9ea16c3f10>
Collect all the outline files into one folder.
The LOCO-EFA method was created by Yara E Sánchez-Corrales, Matthew Hartley, Jop van Rooij, Athanasius F M Marée, Verônica A Grieneisen (Development. 2018 Mar 20;145(6):dev156778. doi: 10.1242/dev.156778.) LOCO-EFA stands for 'Lobe contribution elliptical Fourier analysis (LOCO-EFA)'. This method generates from digitalised two-dimensional cell outlines meaningful descriptors that can be directly matched to morphological features.
Link to paper: https://pubmed.ncbi.nlm.nih.gov/29444894/
Their (example) code can be found here: https://bitbucket.org/mareelab/loco_efa/src/master/src/loco-efa_example.c
We slightly adjusted the code so that the L-values found for each outline were stored in a text file in a path that we set in line 17 in our .c file. When using our file, you should adjust this path to a path on your pc. Name of our file: loco-efa_example_GS.c
Below's code worked on our MacBook Pro from 2017.
First, the code needs to be compiled (C is a mid-level language and it needs a compiler to convert it into an executable code so that the program can be run on our machine).
<br> a. Check the macOS version of your system and check online the release date of this version.
<br> b. Go to https://developer.apple.com/download/all/?q=command%20line%20tools and download the command line tools for Xcode that were released before the release date of your latest macOS version.
c. Install the command line tools
# Note: some of our images were overlapping each other, which caused some gastruloids to appear twice in our data set.
# We carefully checked all the extracted outlines and removed the doubles. Therefore, 500 - 413 = 87 extracted outlines were removed.
# 96 h gastruloids
a1 = np.arange(1,7)
a2 = np.arange(10,27)
a3 = np.arange(29,37)
a4 = np.array([38, 40])
a5 = np.arange(42,90)
a6 = np.arange(91,100)
a62 = np.array([102,104,106,107,109,112,118,121])
a7 = np.arange(123,157)
a8 = np.array([158,166])
a9 = np.arange(168,172)
a10 = np.array([173])
a11 = np.arange(175,183)
a12 = np.array([185,186,188])
a13 = np.arange(192,214)
a14 = np.arange(216,224)
a15 = np.array([225,226])
a16 = np.arange(228,243)
a17 = np.arange(244,252)
a18 = np.array([253,254])
a19 = np.arange(256,273)
a20 = np.array([274])
a21 = np.arange(276,282)
a22 = np.array([283])
a23 = np.arange(285,296)
a24 = np.array([297])
a25 = np.arange(299,304)
a26 = np.arange(306,311)
a27 = np.array([312,317,319,322,323])
a28 = np.arange(326,349)
a29 = np.arange(350,363)
a30 = np.arange(364,372)
a31 = np.arange(373,376)
a32 = np.arange(377,380)
a33 = np.arange(381,390)
a34 = np.arange(391,400)
a35 = np.arange(401,408)
a36 = np.arange(409,412)
a37 = np.arange(413,435)
a38 = np.arange(437,448)
a39 = np.arange(449,452)
a40 = np.arange(453,459)
a41 = np.arange(461,472)
a42 = np.array([473,474])
a43 = np.arange(476,485)
a44 = np.arange(486,490)
a45 = np.arange(492,500)
a_96h = np.concatenate((a1,a2,a3,a4,a5,a6,a62,a7,a8,a9,a10,a11,a12,a13,a14,a15,a16,a17,a18,a19,a20,a21,a22,a23,a24,a25,a26,a27,a28,a29,a30,a31,a32,a33,a34,a35,a36,a37,a38,a39,a40,a41,a42,a43,a44,a45),dtype = np.int64)
print('Total number of 96 h gastruloids:', len(a_96h))
# 72 h gastruloids
a1 = np.arange(1,11)
a2 = np.arange(13,27)
a3 = np.array([30])
a4 = np.arange(32,89)
a5 = np.arange(91,101)
a6 = np.array([112])
a7 = np.array([114])
a8 = np.arange(116,127)
a9 = np.arange(130,134)
a10 = np.arange(136,140)
a11 = np.arange(146,147)
a12 = np.arange(149,155)
a13 = np.arange(159,171)
a_72h = np.concatenate((a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13))#,dtype = np.int64)
print('Total number of 72 h gastruloids:', len(a_72h))
# 10 uM ROCK inhibitor gastruloids
a1 = np.arange(1,19)
a2 = np.array([20])
a3 = np.arange(23,37)
a4 = np.array([38])
a5 = np.array([40])
a6 = np.arange(42,73)
a7 = np.arange(74,76)
a8 = np.arange(77,112)
a9 = np.array([113])
a10 = np.arange(115,118)
a11 = np.array([119])
a12 = np.array([121])
a13 = np.array([123])
a_10uM_RI = np.concatenate((a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13))#,dtype = np.int64)
print('ROCK inhibitor experiment - total number of treated gastruloids:', len(a_10uM_RI))
# No ROCK inhibitor gastruloids
a1 = np.arange(1,30)
a2 = np.arange(31,34)
a3 = np.arange(35,50)
a4 = np.arange(53,58)
a5 = np.arange(60,67)
a6 = np.arange(68,80)
a7 = np.arange(81,96)
a8 = np.arange(97,102)
a9 = np.arange(103,105)
a10 = np.arange(106,113)
a11 = np.array([115])
a12 = np.arange(117,122)
a_no_RI = np.concatenate((a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12))#,dtype = np.int64)
print('ROCK inhibitor experiment - total number of control gastruloids:', len(a_no_RI))
Total number of 96 h gastruloids: 413 Total number of 72 h gastruloids: 132 ROCK inhibitor experiment - total number of treated gastruloids: 110 ROCK inhibitor experiment - total number of control gastruloids: 106
folder_72h = "/Users/esmeeadegeest/Documents/Outlines_72h_143unique/Outlines"
L2_div_L1_72h_array=[]
L3_div_L1_72h_array=[]
for i in range(np.size(a_72h)):
i_string = str(a_72h[i])
data_file = "Contour_"+i_string+"outline.txt"
path_i = os.path.join(folder_72h,data_file)
contents = pd.read_table(path_i)
L_values_72h = pd.DataFrame(contents).to_numpy() # convert dataframe to array in numpy
L2_div_L1_72h = np.divide(L_values_72h[1], L_values_72h[0]) # divide L2 by L1
L3_div_L1_72h = np.divide(L_values_72h[2], L_values_72h[0]) # divide L3 by L1
L2_div_L1_72h_array.append(L2_div_L1_72h) # store L2/L1 values
L3_div_L1_72h_array.append(L3_div_L1_72h) # store L3/L1 values
folder_96h = "/Users/esmeeadegeest/Documents/L_values_no_doubles_413/"
L2_div_L1_96h_array=[]
L3_div_L1_96h_array=[]
for i in range(np.size(a_96h)):
i_string = str(a_96h[i])
data_file = "Data_"+i_string+".txt"
path_i = os.path.join(folder_96h,data_file)
contents = pd.read_table(path_i)
L_values_96h = pd.DataFrame(contents).to_numpy() # convert dataframe to array in numpy
L2_div_L1_96h = np.divide(L_values_96h[1], L_values_96h[0]) # divide L2 by L1
L3_div_L1_96h = np.divide(L_values_96h[2], L_values_96h[0]) # divide L3 by L1
L2_div_L1_96h_array.append(L2_div_L1_96h) # store L2/L1 values
L3_div_L1_96h_array.append(L3_div_L1_96h) # store L3/L1 values
folder_No_RI = "/Users/esmeeadegeest/Documents/L-values_No-RI/"
L2_div_L1_No_RI_array=[]
L3_div_L1_No_RI_array=[]
for i in range(np.size(a_no_RI)):
i_string = str(a_no_RI[i])
data_file = "Contour_"+i_string+"outline.txt"
path_i = os.path.join(folder_No_RI,data_file)
contents = pd.read_table(path_i)
L_values_No_RI = pd.DataFrame(contents).to_numpy() # convert dataframe to array in numpy
L2_div_L1_No_RI = np.divide(L_values_No_RI[1], L_values_No_RI[0]) # divide L2 by L1
L3_div_L1_No_RI = np.divide(L_values_No_RI[2], L_values_No_RI[0]) # divide L3 by L1
L2_div_L1_No_RI_array.append(L2_div_L1_No_RI) # store L2/L1 values
L3_div_L1_No_RI_array.append(L3_div_L1_No_RI) # store L3/L1 values
folder_10uM_RI = "/Users/esmeeadegeest/Documents/L-values_10uM-RI/"
L2_div_L1_10uM_RI_array=[]
L3_div_L1_10uM_RI_array=[]
for i in range(np.size(a_10uM_RI)):
i_string = str(a_10uM_RI[i])
data_file = "Contour_"+i_string+"outline.txt"
path_i = os.path.join(folder_10uM_RI,data_file)
contents = pd.read_table(path_i)
L_values_10uM_RI = pd.DataFrame(contents).to_numpy() # convert dataframe to array in numpy
L2_div_L1_10uM_RI = np.divide(L_values_10uM_RI[1], L_values_10uM_RI[0]) # divide L2 by L1
L3_div_L1_10uM_RI = np.divide(L_values_10uM_RI[2], L_values_10uM_RI[0]) # divide L3 by L1
L2_div_L1_10uM_RI_array.append(L2_div_L1_10uM_RI) # store L2/L1 values
L3_div_L1_10uM_RI_array.append(L3_div_L1_10uM_RI) # store L3/L1 values
#72 h and 96 h
%matplotlib inline
font = {'size' : 12}
matplotlib.rc('font', **font)
matplotlib.rcParams['font.sans-serif'] = "Arial"
matplotlib.rcParams['font.family'] = "sans-serif"
cbmap=['#377eb8', '#ff7f00', '#4daf4a',
'#f781bf', '#a65628', '#984ea3',
'#999999', '#e41a1c', '#dede00']
plt.gcf().clear()
fig = plt.figure(1)
ax = fig.add_subplot(111)
ax.scatter(L2_div_L1_72h_array,L3_div_L1_72h_array, s=1, marker ="o", label='72h gastruloids', c=cbmap[5])
ax.scatter(L2_div_L1_96h_array,L3_div_L1_96h_array, s=1, marker ="o", label='96h gastruloids', c=cbmap[1])
fig.set_size_inches(10/2.54, 8/2.54)
plt.xlabel('$L_2/L_1$')
plt.ylabel('$L_3/L_1$')
#plt.title('all <-> all')
handles, labels = ax.get_legend_handles_labels()
lgd = ax.legend(handles, labels, loc='upper left', bbox_to_anchor=(1,1))
fig.set_size_inches(10/2.54, 8/2.54)
# save figure
#fig.savefig('Figures_paper/Figure_1_LOCO-EFA_72h_96h_purple_orange.svg', bbox_extra_artists=([lgd]), bbox_inches='tight')
# ROCK inhibitor experiment
%matplotlib inline
font = {'size' : 12}
matplotlib.rc('font', **font)
matplotlib.rcParams['font.sans-serif'] = "Arial"
matplotlib.rcParams['font.family'] = "sans-serif"
cbmap=['#377eb8', '#ff7f00', '#4daf4a',
'#f781bf', '#a65628', '#984ea3',
'#999999', '#e41a1c', '#dede00']
plt.gcf().clear()
fig = plt.figure(1)
ax = fig.add_subplot(111)
ax.scatter(L2_div_L1_No_RI_array,L3_div_L1_No_RI_array, s=1, marker ="o", label='Control', c=cbmap[1])
ax.scatter(L2_div_L1_10uM_RI_array,L3_div_L1_10uM_RI_array, s=1, marker ="o", label=r'10 $\mu$M Y-27632', c=cbmap[5])
fig.set_size_inches(10/2.54, 8/2.54)
plt.xlabel('$L_2/L_1$')
plt.ylabel('$L_3/L_1$')
plt.xlim(-0.02, 0.7)
plt.ylim(-0.015, 0.305)
#plt.title('all <-> all')
handles, labels = ax.get_legend_handles_labels()
lgd = ax.legend(handles, labels, loc='upper left', bbox_to_anchor=(1,1))
fig.set_size_inches(10/2.54, 8/2.54)
# save figure
#fig.savefig('Documents/Figure_X_LOCO-EFA_RI_purple_orange.svg', bbox_extra_artists=([lgd]), bbox_inches='tight')
Our time lapse data consists of 31 image stacks, corresponding to 31 time points. Here, we used maximum projections in Z to analyze the 3D data in 2D. To determine where cells are with respect to the long and short axis of the gastruloid, we fit an ellipse to the 2D gastruloid shape. To obtain the 2D gastruloid shape, we first pre-process the images in FIJI to enhance the contrast. Then we use several morphological operations to create a binary image that shows the shape of the gastruloid. The binary image is a rough estimate of the gastruloid's shape, since not all cells within the gastruloid were mCherry+. From the binary images, an outline is created, which is used to fit an ellipse to.
# Pre-processing in FIJI:
# run("RGB Color");
# run("Enhance Contrast", "saturated=0.35");
# should be stored as PNG file.)
%matplotlib inline
img = []
img_GB = []
img_GB_thresh = []
img_binary = []
plt.rcParams['figure.figsize'] = [30, 20]
kernel = np.ones((5, 5), np.uint8) # Taking a matrix of size 5 as the kernel for the morphological operations
# Perform morphological operations on the maximum z projection for each time point using a specified threshold value and store the images in img_binary
for i in range(1,32):
img_read = cv2.imread('Mosaic_TL_mCherry/Max_projections_Z/MAX_T'+str(i)+'_TL-P1-10min-33ts_G2_561nm-1700mA_160ms_6p7_4um_40zimage1_crop.png',0)
img.append(img_read)
if i == 1:
threshold_value = 63 # We used 40 for all frames, except for frame of t = 1 (63) and t = 26 (45).
elif i == 26:
threshold_value = 45
else:
threshold_value = 40
# Apply Gaussian blur
img_read_GB = cv2.GaussianBlur(img_read, (51,51), cv2.BORDER_DEFAULT) # kernel size should be odd
img_GB.append(img_read_GB)
# Apply threshold value >> check in the next cell for each time frame whether the threshold value is correct
ret,img_GB_read_thresh1 = cv2.threshold(img_read_GB,threshold_value,255,cv2.THRESH_BINARY)
img_GB_thresh.append(img_GB_read_thresh1)
# Apply dilation multiple times to fill the holes.
img_dilation = cv2.dilate(img_GB_read_thresh1, kernel, iterations=40)
# Apply subsequently erosion to retrieve the initial size
img_erosion = cv2.erode(img_dilation, kernel, iterations=40)
img_binary.append(img_erosion)
%matplotlib inline
plt.rcParams['figure.figsize'] = [15, 15]
fig, axs = plt.subplots(2, 2)
# Show for one time point the original image, the image after Gaussian blur,
# the image after applying a threshold, and the image after dilation and erosion.
n = 31 # time point
axs[0,0].imshow(img[n-1], cmap='gray')
axs[0,1].imshow(img_GB[n-1], cmap='gray')
axs[1,0].imshow(img_GB_thresh[n-1], cmap='gray')
axs[1,1].imshow(img_binary[n-1], cmap='gray')
plt.show()
# Save binary images, note down the threshold value for each image
# We used threshold value 40 for all frames, except for the frame of t = 1 (63) and t = 26 (45).
# t = 1, threshold = 63
cv2.imwrite('Mosaic_TL_mCherry/Binary_images/t1-63-binary.png', img_binary[0])
# t = 2-25, threshold = 40
for i in range(2,26):
cv2.imwrite('Mosaic_TL_mCherry/Binary_images/t'+str(i)+'-40-binary.png', img_binary[i-1])
# t = 26, threshold = 45
cv2.imwrite('Mosaic_TL_mCherry/Binary_images/t26-45-binary.png', img_binary[25])
# t = 27-31, threshold = 40
for i in range(27,32):
cv2.imwrite('Mosaic_TL_mCherry/Binary_images/t'+str(i)+'-40-binary.png', img_binary[i-1])
# Area calculation and comparison t1 and t31
image_loc_t1 = 'Mosaic_TL_mCherry/Binary_images/t1-63-binary.png'
image_loc_t31 = 'Mosaic_TL_mCherry/Binary_images/t31-40-binary.png'
mask_t1 = cv2.imread(image_loc_t1)
mask_t31 = cv2.imread(image_loc_t31)
imgray_t1 = cv2.cvtColor(mask_t1,cv2.COLOR_BGR2GRAY)
imgray_t31 = cv2.cvtColor(mask_t31,cv2.COLOR_BGR2GRAY)
area_t1 = np.sum(imgray_t1)/255
area_t31 = np.sum(imgray_t31)/255
print('Ratio area t31/t1:', area_t31/area_t1)
plt.rcParams['figure.figsize'] = [10, 10]
fig, axs = plt.subplots(1, 2)
axs[0].imshow(mask_t1, cmap='gray')
axs[1].imshow(mask_t31, cmap='gray')
Ratio area t31/t1: 1.3836912100209595
<matplotlib.image.AxesImage at 0x7fb20addab80>
# We used the following code to extract outlines from the binary images.
def extract_outline_iv():
for j in range(1,32):
ind = list() # index
cont_all = list() # contours
if j == 1:
image_loc = 'Mosaic_TL_mCherry/Binary_images/t'+str(j)+'-63-binary.png'
file_loc = 'Mosaic_TL_mCherry/Binary_outlines/t'+str(j)+'-63-binary_outline.txt'
elif j == 26:
image_loc = 'Mosaic_TL_mCherry/Binary_images/t'+str(j)+'-45-binary.png'
file_loc = 'Mosaic_TL_mCherry/Binary_outlines/t'+str(j)+'-45-binary_outline.txt'
else:
image_loc = 'Mosaic_TL_mCherry/Binary_images/t'+str(j)+'-40-binary.png'
file_loc = 'Mosaic_TL_mCherry/Binary_outlines/t'+str(j)+'-40-binary_outline.txt'
file1 = open(file_loc,"a")
mask = cv2.imread(image_loc)
if mask is not None:
imgray = cv2.cvtColor(mask,cv2.COLOR_BGR2GRAY)
ret,thresh = cv2.threshold(imgray,200,255,0)
contours, im2 = cv2.findContours(thresh,cv2.RETR_TREE,cv2.CHAIN_APPROX_NONE)
maxsizeloc = 0
maxsize = 0
for i in range (len(contours)):
if contours[i].shape[0] > maxsize:
maxsizeloc = i
maxsize = contours[i].shape[0]
cont = np.zeros((contours[maxsizeloc].shape[0], 2))
cont[:,0] = contours[maxsizeloc][:,0,0]
cont[:,1] = contours[maxsizeloc][:,0,1]
plt.plot(cont[:,0], cont[:,1])
cont_all.append(cont)
# Write outlines in txt file
for i in range (maxsize):
file1.write(str(cont_all[0][i][0]) + ',' + str(cont_all[0][i][1]) + '\n')
file1.close()
return(cont_all)
extract_outline_iv()
# Ellipse fit parameters
xc_iv = np.empty(31) # center coordinate x of ellipse
yc_iv = np.empty(31) # center coordinate y of ellipse
a_iv = np.empty(31) # length short arm of ellipse
b_iv = np.empty(31) # length long arm of ellipse
theta_iv= np.empty(31) # rotation angle of ellipse
# Figure
%matplotlib inline
fig, axs = plt.subplots(2, 16,figsize=(100,30))#sharex=True, sharey=True)
# Load outline coordinates for each time step and fit ellipse, store ellipse parameters
for j in range (1,32):
# txt file name in which outline is stored
if j == 1:
outline_loc = 'Mosaic_TL_mCherry/Binary_outlines/t'+str(j)+'-63-binary_outline.txt'
file_loc = 'Mosaic_TL_mCherry/Binary_outlines/t'+str(j)+'-63_xycenter-theta-axes.txt' # assign a file location and open it to later on write the ellipse parameters in the file
elif j == 26:
outline_loc = 'Mosaic_TL_mCherry/Binary_outlines/t'+str(j)+'-45-binary_outline.txt'
file_loc = 'Mosaic_TL_mCherry/Binary_outlines/t'+str(j)+'-45_xycenter-theta-axes.txt'
else:
outline_loc = 'Mosaic_TL_mCherry/Binary_outlines/t'+str(j)+'-40-binary_outline.txt'
file_loc = 'Mosaic_TL_mCherry/Binary_outlines/t'+str(j)+'-40_xycenter-theta-axes.txt'
file1 = open(file_loc,"a")
# open outline file and store coordinates in array
with open(outline_loc) as file_name:
points = np.loadtxt(file_name, delimiter=",")
a_points = np.array(points)
x = a_points[:, 0]
y = a_points[:, 1]
# fit ellipse to array and get parameters
ell = EllipseModel()
ell.estimate(a_points)
xc_iv[j-1], yc_iv[j-1], a_iv[j-1], b_iv[j-1], theta_iv[j-1] = ell.params
# Write parameters to txt file
#file1.write(str(xc_iv[j-1]) + ',' + str(yc_iv[j-1]) + ',' + str(theta_iv[j-1]) + ',' + str(a_iv[j-1]) + ',' + str(b_iv[j-1]) + '\n')
# Plot ellipse fits
if j < 16: # top row of images
axs[0,j-1].scatter(x, y) # outline
axs[0,j-1].scatter(xc_iv[j-1], yc_iv[j-1], color='red', s=100) # center of ellipse
axs[0,j-1].set_xlim(0,1600)
axs[0,j-1].set_ylim(0,3000)
axs[0,j-1].set_xticks([])
axs[0,j-1].set_yticks([])
# ellipse
ell_patch = Ellipse((xc_iv[j-1], yc_iv[j-1]), 2*a_iv[j-1], 2*b_iv[j-1], theta_iv[j-1]*180/np.pi, edgecolor='red', facecolor='none')
# ellipse with orientation angle 0
ell_patch2 = Ellipse((xc_iv[j-1], yc_iv[j-1]), 2*a_iv[j-1], 2*b_iv[j-1], 0, edgecolor='black', facecolor='none')
axs[0,j-1].add_patch(ell_patch)
axs[0,j-1].add_patch(ell_patch2)
else: # bottom row of images
axs[1,j-16].scatter(x, y) # outline
axs[1,j-16].scatter(xc_iv[j-1], yc_iv[j-1], color='red', s=100) # center of ellipse
axs[1,j-16].set_xlim(0,1600)
axs[1,j-16].set_ylim(0,3000)
axs[1,j-16].set_xticks([])
axs[1,j-16].set_yticks([])
# ellipse
ell_patch = Ellipse((xc_iv[j-1], yc_iv[j-1]), 2*a_iv[j-1], 2*b_iv[j-1], theta_iv[j-1]*180/np.pi, edgecolor='red', facecolor='none')
# ellipse with orientation angle 0
ell_patch2 = Ellipse((xc_iv[j-1], yc_iv[j-1]), 2*a_iv[j-1], 2*b_iv[j-1], 0, edgecolor='black', facecolor='none')
axs[1,j-16].add_patch(ell_patch)
axs[1,j-16].add_patch(ell_patch2)
plt.show()
# From the images we can tell that for some shapes the short arm length and long arm length are
# stored in the wrong parameters: b and a, respectively, while they should be stored in a and b, respectively.
# Here we print the short arm length (a) and long arm length (b) and see that they are not always stored correctly.
print('Time point 1 (correct): a = ',a_iv[0],'b =', b_iv[0])
print('Time point 4 (incorrect): a = ',a_iv[3],'b = ', b_iv[3])
Time point 1 (correct): a = 588.7241640888858 b = 879.3560145828596 Time point 4 (incorrect): a = 921.3089122428711 b = 591.6155040634976
# We traced cells manually using the point tool in ImageJ. The trajectories were stored in an ROI set.
# Import the roi_set of each cell
N_cells = 14
cell_list = [None] * N_cells
roi_set = [None] * N_cells
for i in range(N_cells):
path_file = 'Mosaic_TL_mCherry/Cells_by_number/RoiSet_Cell_' + str(i+1) + '.zip'
roi_set[i] = read_roi_zip(path_file)
We also had dividing cells and therefore tracked the daughter cells after division. One cell was tracked starting from the last time point to the first time point for convenience. The following cells have different scripts to convert their ROIs to arrays:
cell 3 = dividing cell, positions of daughter cells alternate
cell 4 = dividing cell, position of first daughter cell first
cell 5 = dividing cell, position of first daughter cell first
cell 7 = dividing cell, position of first daughter cell first
cell 8 = order in set is from last to first time point
cell 9 = dividing cell, position of first daughter cell first
For each cell, we store the trajectories in a list. Depending on whether a cell divides, we store the trajectory of the mother cell in path_x, path_y and the trajectories of the daughter cells in path_x_d1, path_y_d1 and path_x_d2, path_y_d2.
cell_list = path_x, path_y, division_status, time_step_division, path_x_d1, path_y_d1, path_x_d2, path_y_d2
# Dividing cell script for cell 3 (alternating position)
n = 3
roi_set_sel = roi_set[n-1]
time_step_division = 17
t = 1 # start time
division_status = 'Dividing cell'
path_x = np.arange(0) # mother cell x
path_y = np.arange(0) # mother cell y
path_x_d1 = np.arange(0) # daughter cell 1 x
path_x_d2 = np.arange(0) # daughter cell 2 x
path_y_d1 = np.arange(0) # daughter cell 1 y
path_y_d2 = np.arange(0) # daughter cell 2 xy
roi_set_length = len(roi_set_sel)
d1 = bool(True) # daughter cell 1 is set to true
for key, value in roi_set_sel.items():
path_roi = roi_set_sel[key]
if t < time_step_division:
path_x = np.append(path_x, path_roi['x'][0])
path_y = np.append(path_y, path_roi['y'][0])
t = t + 1
elif t >= time_step_division:
if d1 == True:
path_x_d1 = np.append(path_x_d1, path_roi['x'][0])
path_y_d1 = np.append(path_y_d1, path_roi['y'][0])
d1 = bool(False) # set daughter cell 1 to false > store next position in array of 2nd daughter cell
t = t + 1
elif d1 == False:
path_x_d2 = np.append(path_x_d2, path_roi['x'][0])
path_y_d2 = np.append(path_y_d2, path_roi['y'][0])
d1 = bool(True) # set daughter cell 1 to true > store next position in array of 1st daughter cell
t = t + 1
# store the trajectories in a list
cell_list[n-1] = [path_x, path_y, division_status, time_step_division, path_x_d1[:(len(path_x_d1))], path_y_d1[:(len(path_y_d1))], path_x_d2[:(len(path_x_d2))], path_y_d2[:(len(path_y_d2))]]
# Dividing cell script for cell 5 (weird order: first alternating positions, then paths subsequently stored in ROI set)
n = 5
roi_set_sel = roi_set[n-1]
time_step_division_par = 9
time_step_division_sep = 16
t_sep_start = time_step_division_par + 2*(time_step_division_sep - time_step_division_par)
t_sep_d2_start = t_sep_start + (31 - time_step_division_sep + 1)
roi_set_length = len(roi_set_sel)
t = 1
division_status = 'Dividing cell'
path_x = np.arange(0)
path_y = np.arange(0)
path_x_d1 = np.arange(0)
path_x_d2 = np.arange(0)
path_y_d1 = np.arange(0)
path_y_d2 = np.arange(0)
d1 = bool(True)
for key, value in roi_set_sel.items():
path_roi = roi_set_sel[key]
if t < time_step_division_par: # to 8 (incl)
path_x = np.append(path_x, path_roi['x'][0])
path_y = np.append(path_y, path_roi['y'][0])
t = t + 1
elif t >= time_step_division_par and t < t_sep_start: # 9 to 15 (incl)
if d1 == True:
path_x_d1 = np.append(path_x_d1, path_roi['x'][0])
path_y_d1 = np.append(path_y_d1, path_roi['y'][0])
d1 = bool(False)
t = t + 1
elif d1 == False:
path_x_d2 = np.append(path_x_d2, path_roi['x'][0])
path_y_d2 = np.append(path_y_d2, path_roi['y'][0])
d1 = bool(True)
t = t + 1
elif t >= t_sep_start and t < t_sep_d2_start:
path_x_d1 = np.append(path_x_d1, path_roi['x'][0])
path_y_d1 = np.append(path_y_d1, path_roi['y'][0])
t = t + 1
elif t >= t_sep_d2_start:
path_x_d2 = np.append(path_x_d2, path_roi['x'][0])
path_y_d2 = np.append(path_y_d2, path_roi['y'][0])
t = t + 1
cell_list[n-1] = [path_x, path_y, division_status, time_step_division_par, path_x_d1[:(len(path_x_d1))], path_y_d1[:(len(path_y_d1))], path_x_d2[:(len(path_x_d2))], path_y_d2[:(len(path_y_d2))]]
# Dividing cell script for cell 7
n = 7
time_step_division = 14
roi_set_sel = roi_set[n-1]
roi_set_length = len(roi_set_sel)
t = 1
division_status = 'Dividing cell'
path_x = np.arange(0)
path_y = np.arange(0)
path_x_d1 = np.arange(0)
path_x_d2 = np.arange(0)
path_y_d1 = np.arange(0)
path_y_d2 = np.arange(0)
for key, value in roi_set_sel.items():
path_roi = roi_set_sel[key]
if t < time_step_division:
path_x = np.append(path_x, path_roi['x'][0])
path_y = np.append(path_y, path_roi['y'][0])
t = t + 1
elif t >= time_step_division and t <= 31:
path_x_d1 = np.append(path_x_d1, path_roi['x'][0])
path_y_d1 = np.append(path_y_d1, path_roi['y'][0])
t = t + 1
elif t > 31:
path_x_d2 = np.append(path_x_d2, path_roi['x'][0])
path_y_d2 = np.append(path_y_d2, path_roi['y'][0])
t = t + 1
cell_list[n-1] = [path_x, path_y, division_status, time_step_division, path_x_d1[:(len(path_x_d1))], path_y_d1[:(len(path_y_d1))], path_x_d2[:(len(path_x_d2))], path_y_d2[:(len(path_y_d2))]]
# Dividing cell script for cell 9
n = 9
time_step_division = 12
roi_set_sel = roi_set[n-1]
roi_set_length = len(roi_set_sel)
t = 1
division_status = 'Dividing cell'
path_x = np.arange(0)
path_y = np.arange(0)
path_x_d1 = np.arange(0)
path_x_d2 = np.arange(0)
path_y_d1 = np.arange(0)
path_y_d2 = np.arange(0)
for key, value in roi_set_sel.items():
path_roi = roi_set_sel[key]
if t < time_step_division:
path_x = np.append(path_x, path_roi['x'][0])
path_y = np.append(path_y, path_roi['y'][0])
t = t + 1
elif t >= time_step_division and t <= 31:
path_x_d1 = np.append(path_x_d1, path_roi['x'][0])
path_y_d1 = np.append(path_y_d1, path_roi['y'][0])
t = t + 1
elif t > 31:
path_x_d2 = np.append(path_x_d2, path_roi['x'][0])
path_y_d2 = np.append(path_y_d2, path_roi['y'][0])
t = t + 1
cell_list[n-1] = [path_x, path_y, division_status, time_step_division, path_x_d1[:(len(path_x_d1))], path_y_d1[:(len(path_y_d1))], path_x_d2[:(len(path_x_d2))], path_y_d2[:(len(path_y_d2))]]
# Non-dividing cells: 1, 2, 4, 6, 10, 11, 12, 13, 14
# Fill in here:
n = [1, 2, 4, 6, 10, 11, 12, 13, 14]
division_status = 'Non-dividing cell'
for i in n:
path_x = np.arange(0)
path_y = np.arange(0)
roi_set_sel = roi_set[i-1]
for key, value in roi_set_sel.items():
path_roi = roi_set_sel[key]
path_x = np.append(path_x,path_roi['x'][0])
path_y = np.append(path_y,path_roi['y'][0])
cell_list[i-1] = [path_x[:(len(path_x))], path_y[:(len(path_y))], division_status]
# Non-dividing cell - cell 8 (reverse order)
# Fill in here:
n = 8
roi_set_sel = roi_set[n-1]
division_status = 'Non-dividing cell'
path_x = np.arange(0)
path_y = np.arange(0)
for key, value in roi_set_sel.items():
path_roi = roi_set_sel[key]
path_x = np.append(path_x,path_roi['x'][0])
path_y = np.append(path_y,path_roi['y'][0])
path_x_rev = path_x[::-1]
path_y_rev = path_y[::-1]
cell_list[n-1] = [path_x_rev[:(len(path_x_rev))], path_y_rev[:(len(path_y_rev))], division_status]
n_cells = 14
divcells = 4 # 4 dividing cells
#cells_sel = [1, 2, 4, 6, 8, 10, 11, 12, 13, 14] #non-dividing cells only
cells_sel = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14] #all cells
div_cells_no = [3, 5, 7, 9]
t_steps = 31
go_twice = False # is true for dividing cells
count = 0
theta_corr_iv = [None] * t_steps
rotated_cell_centers_iv = np.empty((n_cells+divcells, t_steps, 4)) # we store: x, y, V_x, V_y
for i in range(n_cells):
sel_cell = cells_sel[i]-1
if cell_list[sel_cell][2] == 'Dividing cell':
go_twice = True
x_d1 = []
y_d1 = []
x_d2 = []
y_d2 = []
x_d1 = np.append([cell_list[sel_cell][0]],[cell_list[sel_cell][4]]) # paste trajectory (x coordinate) daughter cell 1 after trajectory mother cell
y_d1 = np.append([cell_list[sel_cell][1]],[cell_list[sel_cell][5]]) # paste trajectory (y coordinate) daughter cell 1 after trajectory mother cell
x_d2 = np.append([cell_list[sel_cell][0]],[cell_list[sel_cell][6]]) # paste trajectory (x coordinate) daughter cell 2 to trajectory mother cell
y_d2 = np.append([cell_list[sel_cell][1]],[cell_list[sel_cell][7]]) # paste trajectory (y coordinate) daughter cell 2 to trajectory mother cell
else:
go_twice = False
if go_twice == False: # non-dividing cell
for j in range(t_steps):
if theta_iv[j] > 1.5708: # if gastruloid (ellipse fit) orientation angle is larger than pi/2, correct the angle by subtracting pi/2
theta_corr_iv[j] = theta_iv[j] - 1.5708
else:
theta_corr_iv[j] = theta_iv[j]
# 1. Calculate the difference in x and y position of a cell w.r.t. the center of the ellipse for each time step
dx = cell_list[sel_cell][0][j]-xc_iv[j]
dy = cell_list[sel_cell][1][j]-yc_iv[j]
# 2. Calculate the difference in x and y position of each cell w.r.t. to the new x and y axes, the axes of the ellipse, where the long axis is
# the new y axis and the short axis is the new x axis. We call them V_x and V_y.
V_x = np.cos(theta_corr_iv[j]) * dx + np.sin(theta_corr_iv[j]) * dy
V_y = np.cos(theta_corr_iv[j]) * dy - np.sin(theta_corr_iv[j]) * dx
# 3. Store V_x and V_y and store the new x and y coordinates: V_x + x(center of ellipse) and V_y + y(center of ellipse)
rotated_cell_centers_iv[i+count,j,0] = V_x + xc_iv[j]
rotated_cell_centers_iv[i+count,j,2] = V_x
rotated_cell_centers_iv[i+count,j,1] = V_y + yc_iv[j]
rotated_cell_centers_iv[i+count,j,3] = V_y
elif go_twice == True: # dividing cell
for j in range(t_steps):
if theta_iv[j] > 1.5708: # if gastruloid (ellipse fit) orientation angle is larger than pi/2, correct the angle by subtracting pi/2
theta_corr_iv[j] = theta_iv[j] - 1.5708
else:
theta_corr_iv[j] = theta_iv[j]
# 1. Calculate the difference in x and y position of a cell w.r.t. the center of the ellipse for each time step
dx = x_d1[j]-xc_iv[j]
dy = y_d1[j]-yc_iv[j]
# 2. Calculate the difference in x and y position of each cell w.r.t. to the new x and y axes, the axes of the ellipse, where the long axis is
# the new y axis and the short axis is the new x axis. We call them V_x and V_y.
V_x = np.cos(theta_corr_iv[j]) * dx + np.sin(theta_corr_iv[j]) * dy
V_y = np.cos(theta_corr_iv[j]) * dy - np.sin(theta_corr_iv[j]) * dx
# 3. Store V_x and V_y and store the new x and y coordinates: V_x + x(center of ellipse) and V_y + y(center of ellipse)
rotated_cell_centers_iv[i+count,j,0] = V_x + xc_iv[j]
rotated_cell_centers_iv[i+count,j,2] = V_x
rotated_cell_centers_iv[i+count,j,1] = V_y + yc_iv[j]
rotated_cell_centers_iv[i+count,j,3] = V_y
count = count + 1 # for 2nd daughter cell
for j in range(t_steps):
if theta_iv[j] > 1.5708: # if gastruloid (ellipse fit) orientation angle is larger than pi/2, correct the angle by subtracting pi/2
theta_corr_iv[j] = theta_iv[j] - 1.5708
else:
theta_corr_iv[j] = theta_iv[j]
# 1. Calculate the difference in x and y position of a cell w.r.t. the center of the ellipse for each time step
dx = x_d2[j]-xc_iv[j]
dy = y_d2[j]-yc_iv[j]
# 2. Calculate the difference in x and y position of each cell w.r.t. to the new x and y axes, the axes of the ellipse, where the long axis is
# the new y axis and the short axis is the new x axis. We call them V_x and V_y.
V_x = np.cos(theta_corr_iv[j]) * dx + np.sin(theta_corr_iv[j]) * dy
V_y = np.cos(theta_corr_iv[j]) * dy - np.sin(theta_corr_iv[j]) * dx
# 3. Store V_x and V_y and store the new x and y coordinates: V_x + x(center of ellipse) and V_y + y(center of ellipse)
rotated_cell_centers_iv[i+count,j,0] = V_x + xc_iv[j]
rotated_cell_centers_iv[i+count,j,2] = V_x
rotated_cell_centers_iv[i+count,j,1] = V_y + yc_iv[j]
rotated_cell_centers_iv[i+count,j,3] = V_y
# T = 1, 14 cells starting points
%matplotlib inline
plt.rcParams['figure.figsize'] = [30, 20]
imgplot= mpimg.imread('Mosaic_TL_mCherry/MAX_prediction-t1-200epochsRGB_SB100um.png')
plt.imshow(imgplot)#, extent=[0, 300, 0, 300])
# Colors
cmap = [#'#000000', #black
'#e69d00', #orange
'#f0e442', #yellow
'#009e73', #bluish green
'#56b4e9', #sky blue
'#cc79a7', #reddish purple
'#d55e00', #vermillion (reddish)
'#0072b2', #blue
#'#ffffff', #white
'#7f7f7f'] #grey
cell_marker = ['$1$','$2$','$3$','$4$','$5$','$6$','$7$','$8$','$9$','$10$','$11$','$12$','$13$','$14$']
n = 14
for i in range(n):
if i >= 9:
multiplying_factor = 1.5 # use a larger icon for cells that start to count from 10
else:
multiplying_factor = 1
i_sel = i
cell_marker_i = cell_marker[i_sel]
if i < 8:
color_i = cmap[i]
if i >= 8 and i < 16:
color_i = cmap[i-8]
# Plot starting points
plt.plot(cell_list[i_sel][0][0],cell_list[i_sel][1][0], color = 'k', marker='o',
linewidth=2, markersize=23) # start point black edge
plt.plot(cell_list[i_sel][0][0],cell_list[i_sel][1][0], color = color_i, marker='o',
linewidth=2, markersize=21) # start point colored
plt.plot(cell_list[i_sel][0][0],cell_list[i_sel][1][0], color = 'w', marker='o',
linewidth=2, markersize=16) # start point white
plt.plot(cell_list[i_sel][0][0],cell_list[i_sel][1][0], color = 'k', marker=cell_marker_i,
linewidth=2, markersize=10*multiplying_factor) # start point number
#get current axes
ax = plt.gca()
#hide x-axis
ax.get_xaxis().set_visible(False)
#hide y-axis
ax.get_yaxis().set_visible(False)
# Save figure
#plt.savefig("Mosaic_TL_mCherry/Figure_in_vitro_tracking_mosaic_icons_t1_more_white.png", dpi=600, bbox_inches='tight',pad_inches = 0)
# Plot figure
fig = plt.show()
# T = 31, 14 cells trajectories
%matplotlib inline
plt.rcParams['figure.figsize'] = [30, 20]
imgplot= mpimg.imread('Mosaic_TL_mCherry/MAX_prediction-t31-200epochsRGB-SB100um.png')
plt.imshow(imgplot)#, extent=[0, 300, 0, 300])
# Colors
cmap = [#'#000000', #black
'#e69d00', #orange
'#f0e442', #yellow
'#009e73', #bluish green
'#56b4e9', #sky blue
'#cc79a7', #reddish purple
'#d55e00', #vermillion (reddish)
'#0072b2', #blue
#'#ffffff', #white
'#7f7f7f'] #grey
cell_marker = ['$1$','$2$','$3$','$4$','$5$','$6$','$7$','$8$','$9$','$10$','$11$','$12$','$13$','$14$']
n = 14
for i in range(n):
if i >= 9:
multiplying_factor = 1.5 # use a larger icon for cells that start to count from 10
else:
multiplying_factor = 1
i_sel = i
cell_marker_i = cell_marker[i_sel]
if i < 8:
color_i = cmap[i]
if i >= 8 and i < 16:
color_i = cmap[i-8]
# dividing cells
# division points = cell 3: 17, cell 5: 16, cell 7: 14, cell 9: 12
if cell_list[i_sel][2] == 'Dividing cell':
cell_division_point = cell_list[i_sel][3]
plt.plot(cell_list[i_sel][4], cell_list[i_sel][5], color = 'k', #linestyle = 'dotted',
linewidth = 3, markersize = 3) # black trajectory "daughter cell 1"
plt.plot(cell_list[i_sel][4], cell_list[i_sel][5], color = color_i,
linewidth=2.5, markersize=2.5) # colored trajectory "daughter cell 1"
plt.plot(cell_list[i_sel][6], cell_list[i_sel][7], color = 'k', #linestyle = 'dotted',
linewidth = 3, markersize = 3) # black trajectory "daughter cell 2"
plt.plot(cell_list[i_sel][6], cell_list[i_sel][7], color = color_i,
linewidth=2.5, markersize=2.5) # colored trajectory "daughter cell 2"
# Trajectory from start to end (mother cell or non-dividing cell)
plt.plot(cell_list[i_sel][0], cell_list[i_sel][1], color = 'k', #linestyle = 'dotted',
linewidth = 3, markersize = 3) # line connecting start and end point black
plt.plot(cell_list[i_sel][0], cell_list[i_sel][1], color = color_i,
linewidth=2.5, markersize=2.5) # line connecting start and end point color
# start points
if cell_list[i_sel][2] == 'Dividing cell':
# x and y coordinates for connecting line mother cell and daughter cell
mother_x = cell_list[i_sel][0][-1]
mother_y = cell_list[i_sel][1][-1]
daughter_1_x = cell_list[i_sel][4][0]
daughter_1_y = cell_list[i_sel][5][0]
daughter_2_x = cell_list[i_sel][6][0]
daughter_2_y = cell_list[i_sel][7][0]
x_points_d1 = [mother_x, daughter_1_x]
x_points_d2 = [mother_x, daughter_2_x]
y_points_d1 = [mother_y, daughter_1_y]
y_points_d2 = [mother_y, daughter_2_y]
# plot connecting line for daughter cell 1
plt.plot(x_points_d1, y_points_d1, color = 'k',# linestyle = 'dotted',
linewidth = 3, markersize = 3) # line connecting start and end point black
plt.plot(x_points_d1, y_points_d1, color = color_i,
linewidth=2.5, markersize=2.5) # line connecting start and end point color
# plot connecting line for daughter cell 2
plt.plot(x_points_d2, y_points_d2, color = 'k', #linestyle = 'dotted',
linewidth = 3, markersize = 3) # line connecting start and end point black
plt.plot(x_points_d2, y_points_d2, color = color_i,
linewidth=2.5, markersize=2.5) # line connecting start and end point color
# start point daughter cell 1
plt.plot(cell_list[i_sel][4][0],cell_list[i_sel][5][0], color = 'k',
marker='o', linewidth=20, markersize=23) # end point black edge
plt.plot(cell_list[i_sel][4][0],cell_list[i_sel][5][0], color = color_i,
marker='o', linewidth=20, markersize=21) # end point colored
plt.plot(cell_list[i_sel][4][0],cell_list[i_sel][5][0], color = 'k',
marker='o', linewidth=20, markersize=16) #14 end point black
plt.plot(cell_list[i_sel][4][0],cell_list[i_sel][5][0], color = 'w', marker=cell_marker_i,
linewidth=10, markersize=10*multiplying_factor) # end point number in white
# start point daughter cell 2
plt.plot(cell_list[i_sel][6][0],cell_list[i_sel][7][0], color = 'k',
marker='o', linewidth=20, markersize=23) # end point black edge
plt.plot(cell_list[i_sel][6][0],cell_list[i_sel][7][0], color = color_i,
marker='o', linewidth=20, markersize=21) # end point colored
plt.plot(cell_list[i_sel][6][0],cell_list[i_sel][7][0], color = 'k',
marker='o', linewidth=20, markersize=16) # end point black
plt.plot(cell_list[i_sel][6][0],cell_list[i_sel][7][0], color = 'w', marker=cell_marker_i,
linewidth=20, markersize=10*multiplying_factor) # end point number
# start points non-dividing cells and mother cells
plt.plot(cell_list[i_sel][0][0],cell_list[i_sel][1][0], color = 'k', marker='o',
linewidth=2, markersize=23) # start point black edge
plt.plot(cell_list[i_sel][0][0],cell_list[i_sel][1][0], color = color_i, marker='o',
linewidth=2, markersize=21) # start point colored
plt.plot(cell_list[i_sel][0][0],cell_list[i_sel][1][0], color = 'w', marker='o',
linewidth=2, markersize=16) # start point white
plt.plot(cell_list[i_sel][0][0],cell_list[i_sel][1][0], color = 'k', marker=cell_marker_i,
linewidth=2, markersize=10*multiplying_factor) # start point number
# end points
end_points = 1
if end_points == 1:
if cell_list[i_sel][2] == 'Dividing cell':
# daughter cell 1
plt.plot(cell_list[i_sel][4][-1],cell_list[i_sel][5][-1], color = 'k',
marker='D', linewidth=20, markersize=22) # end point black edge
plt.plot(cell_list[i_sel][4][-1],cell_list[i_sel][5][-1], color = color_i,
marker='D', linewidth=20, markersize=20) # end point colored
plt.plot(cell_list[i_sel][4][-1],cell_list[i_sel][5][-1], color = 'w',
marker='D', linewidth=20, markersize=15) # end point white
plt.plot(cell_list[i_sel][4][-1],cell_list[i_sel][5][-1], color = 'k', marker=cell_marker_i,
linewidth=20, markersize=8*multiplying_factor) # end point number
# daughter cell 2
plt.plot(cell_list[i_sel][6][-1],cell_list[i_sel][7][-1], color = 'k',
marker='D', linewidth=20, markersize=22) # end point black edge
plt.plot(cell_list[i_sel][6][-1],cell_list[i_sel][7][-1], color = color_i,
marker='D', linewidth=20, markersize=20) # end point colored
plt.plot(cell_list[i_sel][6][-1],cell_list[i_sel][7][-1], color = 'w',
marker='D', linewidth=20, markersize=15) #13 end point white
plt.plot(cell_list[i_sel][6][-1],cell_list[i_sel][7][-1], color = 'k', marker=cell_marker_i,
linewidth=20, markersize=10*multiplying_factor) # end point number
else: # no dividing cell
plt.plot(cell_list[i_sel][0][-1],cell_list[i_sel][1][-1], color = 'k',
marker='D', linewidth=20, markersize=22) # end point black edge
plt.plot(cell_list[i_sel][0][-1],cell_list[i_sel][1][-1], color = color_i,
marker='D', linewidth=20, markersize=20) # end point colored
plt.plot(cell_list[i_sel][0][-1],cell_list[i_sel][1][-1], color = 'w',
marker='D', linewidth=20, markersize=15) # end point white
plt.plot(cell_list[i_sel][0][-1],cell_list[i_sel][1][-1], color = 'k', marker=cell_marker_i,
linewidth=20, markersize=10*multiplying_factor) # end point number
#get current axes
ax = plt.gca()
#hide x-axis
ax.get_xaxis().set_visible(False)
#hide y-axis
ax.get_yaxis().set_visible(False)
# Save figure
#plt.savefig("Mosaic_TL_mCherry/Figure_in_vitro_tracking_mosaic_icons-t31_continuous-line_black.png", dpi=600, bbox_inches='tight',pad_inches = 0)
# Plot figure
fig = plt.show()
# CORRECTION FOR DRIFT only
# For T = 31, plot 14 cell trajectories that are corrected for the 'drift' of the center's gastruloid (center of ellipse) over time
# With drift we mean: how the center of the gastruloid - estimated by the center of the fitted ellipse - changes in position over time w.r.t. to the center of the gastruloid at the final time step.
%matplotlib inline
plt.rcParams['figure.figsize'] = [30, 20]
imgplot= mpimg.imread('Mosaic_TL_mCherry/MAX_prediction-t31-200epochsRGB-SB100um.png')
plt.imshow(imgplot)#, extent=[0, 300, 0, 300])
cmap = [#'#000000', #black
'#e69d00', #orange
'#f0e442', #yellow
'#009e73', #bluish green
'#56b4e9', #sky blue
'#cc79a7', #reddish purple
'#d55e00', #vermillion (reddish)
'#0072b2', #blue
#'#ffffff', #white
'#7f7f7f'] #grey
cell_marker = ['$1$','$2$','$3$','$4$','$5$','$6$','$7$','$8$','$9$','$10$','$11$','$12$','$13$','$14$']
n = 14
# Calculate drift
x_corr_drift_iv = [None] * t_steps
y_corr_drift_iv = [None] * t_steps
for j in range(0,t_steps):
x_corr_drift_iv[j] = xc_iv[-1]-xc_iv[j]
y_corr_drift_iv[j] = yc_iv[-1]-yc_iv[j]
# Plot trajectories
for i in range(n):
if i >= 9:
multiplying_factor = 1.5 # use a larger icon for cells that start to count from 10
else:
multiplying_factor = 1
i_sel = i
cell_marker_i = cell_marker[i_sel]
if i < 8:
color_i = cmap[i]
if i >= 8 and i < 16:
color_i = cmap[i-8]
#Trajectories
# Dividing cells
# division points; cell 3: 17, cell 5: 16, cell 7: 14, cell 9: 12
if cell_list[i_sel][2] == 'Dividing cell':
cell_division_point = cell_list[i_sel][3]
plt.plot(cell_list[i_sel][4]+x_corr_drift_iv[-len(cell_list[i_sel][4]):], cell_list[i_sel][5]+y_corr_drift_iv[-len(cell_list[i_sel][5]):], color = 'k',
linewidth = 3, markersize = 3) # black trajectory "daughter cell 1"
plt.plot(cell_list[i_sel][4]+x_corr_drift_iv[-len(cell_list[i_sel][4]):], cell_list[i_sel][5]+y_corr_drift_iv[-len(cell_list[i_sel][5]):], color = color_i,
linewidth=2.5, markersize=2.5) # colored trajectory "daughter cell 1"
plt.plot(cell_list[i_sel][6]+x_corr_drift_iv[-len(cell_list[i_sel][6]):], cell_list[i_sel][7]+y_corr_drift_iv[-len(cell_list[i_sel][7]):], color = 'k',
linewidth = 3, markersize = 3) # black trajectory "daughter cell 2"
plt.plot(cell_list[i_sel][6]+x_corr_drift_iv[-len(cell_list[i_sel][6]):], cell_list[i_sel][7]+y_corr_drift_iv[-len(cell_list[i_sel][7]):], color = color_i,
linewidth=2.5, markersize=2.5) # colored trajectory "daughter cell 2"
# Mother cell
# Trajectory from start to end
plt.plot(cell_list[i_sel][0]+x_corr_drift_iv[0:len(cell_list[i_sel][0])], cell_list[i_sel][1]+y_corr_drift_iv[0:len(cell_list[i_sel][1])], color = 'k',
linewidth = 3, markersize = 3) # line connecting start and end point black
plt.plot(cell_list[i_sel][0]+x_corr_drift_iv[0:len(cell_list[i_sel][0])], cell_list[i_sel][1]+y_corr_drift_iv[0:len(cell_list[i_sel][1])], color = color_i,
linewidth=2.5, markersize=2.5) # line connecting start and end point
# Dividing cells - connecting line and starting points
if cell_list[i_sel][2] == 'Dividing cell':
# Connecting line mother cell and daughter cell
mother_x = cell_list[i_sel][0][-1]+x_corr_drift_iv[len(cell_list[i_sel][0])-1]
mother_y = cell_list[i_sel][1][-1]+y_corr_drift_iv[len(cell_list[i_sel][1])-1]
daughter_1_x = cell_list[i_sel][4][0]+x_corr_drift_iv[len(cell_list[i_sel][0])]
daughter_1_y = cell_list[i_sel][5][0]+y_corr_drift_iv[len(cell_list[i_sel][1])]
daughter_2_x = cell_list[i_sel][6][0]+x_corr_drift_iv[len(cell_list[i_sel][0])]
daughter_2_y = cell_list[i_sel][7][0]+y_corr_drift_iv[len(cell_list[i_sel][1])]
x_points_d1 = [mother_x, daughter_1_x]
x_points_d2 = [mother_x, daughter_2_x]
y_points_d1 = [mother_y, daughter_1_y]
y_points_d2 = [mother_y, daughter_2_y]
plt.plot(x_points_d1, y_points_d1, color = 'k',# linestyle = 'dotted',
linewidth = 3, markersize = 3) # line connecting start and end point black
plt.plot(x_points_d1, y_points_d1, color = color_i,
linewidth=2.5, markersize=2.5) # line connecting start and end point
plt.plot(x_points_d2, y_points_d2, color = 'k', #linestyle = 'dotted',
linewidth = 3, markersize = 3) # line connecting start and end point black
plt.plot(x_points_d2, y_points_d2, color = color_i,
linewidth=2.5, markersize=2.5) # line connecting start and end point
# Start point daughter cell 1
plt.plot(cell_list[i_sel][4][0]+x_corr_drift_iv[len(cell_list[i_sel][0])],cell_list[i_sel][5][0]+y_corr_drift_iv[len(cell_list[i_sel][1])], color = 'k',
marker='o', linewidth=20, markersize=23) # end point black edge
plt.plot(cell_list[i_sel][4][0]+x_corr_drift_iv[len(cell_list[i_sel][0])],cell_list[i_sel][5][0]+y_corr_drift_iv[len(cell_list[i_sel][1])], color = color_i,
marker='o', linewidth=20, markersize=21) # end point colored
plt.plot(cell_list[i_sel][4][0]+x_corr_drift_iv[len(cell_list[i_sel][0])],cell_list[i_sel][5][0]+y_corr_drift_iv[len(cell_list[i_sel][1])], color = 'k',
marker='o', linewidth=20, markersize=16) #14 end point black
plt.plot(cell_list[i_sel][4][0]+x_corr_drift_iv[len(cell_list[i_sel][0])],cell_list[i_sel][5][0]+y_corr_drift_iv[len(cell_list[i_sel][1])], color = 'w', marker=cell_marker_i,
linewidth=10, markersize=10*multiplying_factor) # end point number in white
# Start point daughter cell 2
plt.plot(cell_list[i_sel][6][0]+x_corr_drift_iv[len(cell_list[i_sel][0])],cell_list[i_sel][7][0]+y_corr_drift_iv[len(cell_list[i_sel][1])], color = 'k',
marker='o', linewidth=20, markersize=23) # end point black edge
plt.plot(cell_list[i_sel][6][0]+x_corr_drift_iv[len(cell_list[i_sel][0])],cell_list[i_sel][7][0]+y_corr_drift_iv[len(cell_list[i_sel][1])], color = color_i,
marker='o', linewidth=20, markersize=21) # end point colored
plt.plot(cell_list[i_sel][6][0]+x_corr_drift_iv[len(cell_list[i_sel][0])],cell_list[i_sel][7][0]+y_corr_drift_iv[len(cell_list[i_sel][1])], color = 'k',
marker='o', linewidth=20, markersize=16) # end point black
plt.plot(cell_list[i_sel][6][0]+x_corr_drift_iv[len(cell_list[i_sel][0])],cell_list[i_sel][7][0]+y_corr_drift_iv[len(cell_list[i_sel][1])], color = 'w', marker=cell_marker_i,
linewidth=20, markersize=10*multiplying_factor) # end point number
# Starting points non-dividing cells
plt.plot(cell_list[i_sel][0][0]+x_corr_drift_iv[0],cell_list[i_sel][1][0]+y_corr_drift_iv[0], color = 'k', marker='o',
linewidth=2, markersize=23) # start point black edge
plt.plot(cell_list[i_sel][0][0]+x_corr_drift_iv[0],cell_list[i_sel][1][0]+y_corr_drift_iv[0], color = color_i, marker='o',
linewidth=2, markersize=21) # start point colored
plt.plot(cell_list[i_sel][0][0]+x_corr_drift_iv[0],cell_list[i_sel][1][0]+y_corr_drift_iv[0], color = 'w', marker='o',
linewidth=2, markersize=16) # start point white
plt.plot(cell_list[i_sel][0][0]+x_corr_drift_iv[0],cell_list[i_sel][1][0]+y_corr_drift_iv[0], color = 'k', marker=cell_marker_i,
linewidth=2, markersize=10*multiplying_factor) # start point number
# End points (all cells)
end_points = 1
if end_points == 1:
if cell_list[i_sel][2] == 'Dividing cell':
# daughter cell 1
plt.plot(cell_list[i_sel][4][-1],cell_list[i_sel][5][-1], color = 'k',
marker='D', linewidth=20, markersize=22) # end point black edge
plt.plot(cell_list[i_sel][4][-1],cell_list[i_sel][5][-1], color = color_i,
marker='D', linewidth=20, markersize=20) # end point colored
plt.plot(cell_list[i_sel][4][-1],cell_list[i_sel][5][-1], color = 'w',
marker='D', linewidth=20, markersize=15) # end point white
plt.plot(cell_list[i_sel][4][-1],cell_list[i_sel][5][-1], color = 'k', marker=cell_marker_i,
linewidth=20, markersize=8*multiplying_factor) # end point number
# daughter cell 2
plt.plot(cell_list[i_sel][6][-1],cell_list[i_sel][7][-1], color = 'k',
marker='D', linewidth=20, markersize=22) # end point black edge
plt.plot(cell_list[i_sel][6][-1],cell_list[i_sel][7][-1], color = color_i,
marker='D', linewidth=20, markersize=20) # end point colored
plt.plot(cell_list[i_sel][6][-1],cell_list[i_sel][7][-1], color = 'w',
marker='D', linewidth=20, markersize=15) #13 end point white
plt.plot(cell_list[i_sel][6][-1],cell_list[i_sel][7][-1], color = 'k', marker=cell_marker_i,
linewidth=20, markersize=10*multiplying_factor) # end point number
else: # no dividing cell
plt.plot(cell_list[i_sel][0][-1],cell_list[i_sel][1][-1], color = 'k',
marker='D', linewidth=20, markersize=22) # end point black edge
plt.plot(cell_list[i_sel][0][-1],cell_list[i_sel][1][-1], color = color_i,
marker='D', linewidth=20, markersize=20) # end point colored
plt.plot(cell_list[i_sel][0][-1],cell_list[i_sel][1][-1], color = 'w',
marker='D', linewidth=20, markersize=15) # end point white
plt.plot(cell_list[i_sel][0][-1],cell_list[i_sel][1][-1], color = 'k', marker=cell_marker_i,
linewidth=20, markersize=10*multiplying_factor) # end point number
#get current axes
ax = plt.gca()
#hide x-axis
ax.get_xaxis().set_visible(False)
#hide y-axis
ax.get_yaxis().set_visible(False)
# Save figure
#plt.savefig("Mosaic_TL_mCherry/Figure_in_vitro_tracking_mosaic_icons-t31_continuous-line_black.png", dpi=600, bbox_inches='tight',pad_inches = 0)
# Plot figure
fig = plt.show()
# CORRECTION FOR DRIFT AND ROTATION
# For T = 31, plot 14 cells trajectories that are corrected for the changing rotation of the gastruloid and the 'drift' of the center's gastruloid (center of ellipse) over time.
# With drift we mean: how the center of the gastruloid - estimated by the center of the fitted ellipse - changes in position over time w.r.t. to the center of the gastruloid at the final time step.
# With rotation we mean: how the rotation of the gastruloid changes w.r.t. the gastruloid's orientation at the final time step. To correct for the rotation we use the rotated cell centers and
# rotate them again to the orientation of the gastruloid at the final time step.
%matplotlib inline
plt.rcParams['figure.figsize'] = [30, 20]
imgplot= mpimg.imread('Mosaic_TL_mCherry/MAX_prediction-t31-200epochsRGB-SB100um.png')
plt.imshow(imgplot)#, extent=[0, 300, 0, 300])
# colors per cell
cmap = [#'#000000', #black
'#e69d00', #orange # 1
'#f0e442', #yellow # 2
'#009e73', #bluish green # 3a
'#009e73', #bluish green # 3b
'#56b4e9', #sky blue # 4
'#cc79a7', #reddish purple # 5a
'#cc79a7', #reddish purple # 5b
'#d55e00', #vermillion (reddish) # 6
'#0072b2', #blue # 7a
'#0072b2', #blue # 7b
#'#ffffff', #white
'#7f7f7f', #grey # 8
'#e69d00', #orange # 9a
'#e69d00', #orange # 9b
'#f0e442', #yellow # 10
'#009e73', #bluish green # 11
'#56b4e9', #sky blue # 12
'#cc79a7', #reddish purple # 13
'#d55e00'] #vermillion (reddish) # 14
cell_marker = ['$1$','$2$','$3$','$3$','$4$','$5$','$5$','$6$','$7$','$7$','$8$','$9$','$9$','$10$','$11$','$12$','$13$','$14$']
# Rotate the rotated cell centers to the orientation of the gastruloid at the final time step
no_cells = len(rotated_cell_centers_iv)
rotated_cell_centers_rot_iv = np.empty((no_cells,t_steps,4))
for i in range(no_cells): # cells
for j in range(t_steps): # time
# O = position of cell
O_x_rot_iv = np.cos(2*np.pi - theta_corr_iv[-1]) * (rotated_cell_centers_iv[i,j,0]-xc_iv[j]) + np.sin(2*np.pi - theta_corr_iv[-1]) * (rotated_cell_centers_iv[i,j,1]-yc_iv[j])
rotated_cell_centers_rot_iv[i,j,0] = O_x_rot_iv + xc_iv[j] # rotated x-coordinate of cell
rotated_cell_centers_rot_iv[i,j,2] = O_x_rot_iv # x component of vector from center of ellipse to O
O_y_rot_iv = np.cos(2*np.pi - theta_corr_iv[-1]) * (rotated_cell_centers_iv[i,j,1]-yc_iv[j]) - np.sin(2*np.pi - theta_corr_iv[-1]) * (rotated_cell_centers_iv[i,j,0]-xc_iv[j])
rotated_cell_centers_rot_iv[i,j,1] = O_y_rot_iv + yc_iv[j] # rotated y-coordinate of cell
rotated_cell_centers_rot_iv[i,j,3] = O_y_rot_iv # y component of vector from center of ellipse to O
# Calculate drift
x_corr_drift_iv = [None] * t_steps
y_corr_drift_iv = [None] * t_steps
for j in range(0,t_steps):
x_corr_drift_iv[j] = xc_iv[-1]-xc_iv[j]
y_corr_drift_iv[j] = yc_iv[-1]-yc_iv[j]
# Plot trajectories corrected for DRIFT (add drift) + corrected for ROTATION (plot rotated-rotated cell centers)
for i in range(no_cells):
# icon number size and color
if i > 12:
multiplying_factor = 1.5 # use a larger icon for cells that start to count from 10
else:
multiplying_factor = 1
cell_marker_i = cell_marker[i]
color_i = cmap[i]
# Plot trajectories (rotated-rotated) and add drift
plt.plot(rotated_cell_centers_rot_iv[i,:,0]+x_corr_drift_iv, rotated_cell_centers_rot_iv[i,:,1]+y_corr_drift_iv, color = 'k', #linestyle = 'dotted',
linewidth = 3, markersize = 3) # line connecting start and end point black
plt.plot(rotated_cell_centers_rot_iv[i,:,0]+x_corr_drift_iv, rotated_cell_centers_rot_iv[i,:,1]+y_corr_drift_iv, color = color_i,
linewidth=2.5, markersize=2.5) # line connecting start and end point
# Plot starting points
for i in range(no_cells):
# icon number size and color
if i > 12:
multiplying_factor = 1.5 # use a larger icon for cells that start to count from 10
else:
multiplying_factor = 1
cell_marker_i = cell_marker[i]
color_i = cmap[i]
# For all cells, plot start icon
plt.plot(rotated_cell_centers_rot_iv[i,0,0]+x_corr_drift_iv[0],rotated_cell_centers_rot_iv[i,0,1]+y_corr_drift_iv[0], color = 'k',
marker='o', linewidth=20, markersize=23) # end point black edge
plt.plot(rotated_cell_centers_rot_iv[i,0,0]+x_corr_drift_iv[0],rotated_cell_centers_rot_iv[i,0,1]+y_corr_drift_iv[0], color = color_i,
marker='o', linewidth=20, markersize=21) # end point colored
plt.plot(rotated_cell_centers_rot_iv[i,0,0]+x_corr_drift_iv[0],rotated_cell_centers_rot_iv[i,0,1]+y_corr_drift_iv[0], color = 'w',
marker='o', linewidth=20, markersize=16) #14 point white
plt.plot(rotated_cell_centers_rot_iv[i,0,0]+x_corr_drift_iv[0],rotated_cell_centers_rot_iv[i,0,1]+y_corr_drift_iv[0], color = 'k', marker=cell_marker_i,
linewidth=10, markersize=10*multiplying_factor) # number in black
# For dividing cells, only plot dividing cell (end point mother cell)
if i == 2: # cell 3
sel_cell = i
div_no = cell_list[sel_cell][3] - 2# 17 > 16
if i == 5: # cell 5
sel_cell = i - 1
div_no = cell_list[sel_cell][3] - 2
if i == 8: # cell 7
sel_cell = i - 2
div_no = cell_list[sel_cell][3] - 2
if i == 11: # cell 9
sel_cell = i - 3
div_no = cell_list[sel_cell][3] - 2
if i == 2 or i == 5 or i == 8 or i == 11:
#End point mother cell
plt.plot(rotated_cell_centers_rot_iv[i,div_no,0]+x_corr_drift_iv[div_no],rotated_cell_centers_rot_iv[i,div_no,1]+y_corr_drift_iv[div_no], color = 'k',
marker='o', linewidth=20, markersize=23) # end point black edge
plt.plot(rotated_cell_centers_rot_iv[i,div_no,0]+x_corr_drift_iv[div_no],rotated_cell_centers_rot_iv[i,div_no,1]+y_corr_drift_iv[div_no], color = color_i,
marker='o', linewidth=20, markersize=21) # end point colored
plt.plot(rotated_cell_centers_rot_iv[i,div_no,0]+x_corr_drift_iv[div_no],rotated_cell_centers_rot_iv[i,div_no,1]+y_corr_drift_iv[div_no], color = 'k',
marker='o', linewidth=20, markersize=16) #14 end point black
plt.plot(rotated_cell_centers_rot_iv[i,div_no,0]+x_corr_drift_iv[div_no],rotated_cell_centers_rot_iv[i,div_no,1]+y_corr_drift_iv[div_no], color = 'w', marker=cell_marker_i,
linewidth=10, markersize=10*multiplying_factor) # end point number in white
# For all cells, plot end icon
plt.plot(rotated_cell_centers_rot_iv[i,-1,0]+x_corr_drift_iv[-1],rotated_cell_centers_rot_iv[i,-1,1]+y_corr_drift_iv[-1], color = 'k',
marker='D', linewidth=20, markersize=23) # end point black edge
plt.plot(rotated_cell_centers_rot_iv[i,-1,0]+x_corr_drift_iv[-1],rotated_cell_centers_rot_iv[i,-1,1]+y_corr_drift_iv[-1], color = color_i,
marker='D', linewidth=20, markersize=21) # end point colored
plt.plot(rotated_cell_centers_rot_iv[i,-1,0]+x_corr_drift_iv[-1],rotated_cell_centers_rot_iv[i,-1,1]+y_corr_drift_iv[-1], color = 'w',
marker='D', linewidth=20, markersize=16) #14 point white
plt.plot(rotated_cell_centers_rot_iv[i,-1,0]+x_corr_drift_iv[-1],rotated_cell_centers_rot_iv[i,-1,1]+y_corr_drift_iv[-1], color = 'k', marker=cell_marker_i,
linewidth=10, markersize=10*multiplying_factor) # number in black
#get current axes
ax = plt.gca()
#hide x-axis
ax.get_xaxis().set_visible(False)
#hide y-axis
ax.get_yaxis().set_visible(False)
# Save figure
#plt.savefig("Mosaic_TL_mCherry/Figure_in_vitro_tracking_mosaic_t31_corrected_drift_rotation.png", dpi=600, bbox_inches='tight',pad_inches = 0)
# Plot figure
fig = plt.show()
# T = 1, 14 cells starting points + yellow connecting lines
plt.rcParams['figure.figsize'] = [30, 20]
imgplot= mpimg.imread('Mosaic_TL_mCherry/MAX_prediction-t1-200epochsRGB_SB100um.png')
plt.imshow(imgplot)#, extent=[0, 300, 0, 300])
cmap = [#'#000000', #black
'#e69d00', #orange
'#f0e442', #yellow
'#009e73', #bluish green
'#56b4e9', #sky blue
'#cc79a7', #reddish purple
'#d55e00', #vermillion (reddish)
'#0072b2', #blue
#'#ffffff', #white
'#7f7f7f'] #grey
cell_marker = ['$1$','$2$','$3$','$4$','$5$','$6$','$7$','$8$','$9$','$10$','$11$','$12$','$13$','$14$']
n = 14
# extra lines
#border
x_points_lines = [cell_list[6][0][0],cell_list[2][0][0],cell_list[1][0][0],cell_list[0][0][0],
cell_list[5][0][0],cell_list[11][0][0],cell_list[13][0][0],cell_list[12][0][0],cell_list[9][0][0],cell_list[6][4][0]]
y_points_lines = [cell_list[6][1][0],cell_list[2][1][0],cell_list[1][1][0],cell_list[0][1][0],
cell_list[5][1][0],cell_list[11][1][0],cell_list[13][1][0],cell_list[12][1][0],cell_list[9][1][0],cell_list[6][1][0]]
#in between
x_points_lines_mid_1 = [cell_list[2][0][0],cell_list[3][0][0],cell_list[5][0][0]]
y_points_lines_mid_1 = [cell_list[2][1][0],cell_list[3][1][0],cell_list[5][1][0]]
x_points_lines_mid_2 = [cell_list[6][0][0],cell_list[7][0][0],cell_list[5][0][0]]
y_points_lines_mid_2 = [cell_list[6][1][0],cell_list[7][1][0],cell_list[5][1][0]]
x_points_lines_mid_3 = [cell_list[9][0][0],cell_list[10][0][0],cell_list[11][0][0]]
y_points_lines_mid_3 = [cell_list[9][1][0],cell_list[10][1][0],cell_list[11][1][0]]
plt.plot(x_points_lines, y_points_lines,'y',linewidth = 5)
plt.plot(x_points_lines_mid_1, y_points_lines_mid_1,'y',linewidth = 5)
plt.plot(x_points_lines_mid_2, y_points_lines_mid_2,'y',linewidth = 5)
plt.plot(x_points_lines_mid_3, y_points_lines_mid_3,'y',linewidth = 5)
for i in range(n):
if i >= 9:
multiplying_factor = 1.5 # use a larger icon for cells that start to count from 10
else:
multiplying_factor = 1
i_sel = i
cell_marker_i = cell_marker[i_sel]
if i < 8:
color_i = cmap[i]
if i >= 8 and i < 16:
color_i = cmap[i-8]
# start points
plt.plot(cell_list[i_sel][0][0],cell_list[i_sel][1][0], color = 'k', marker='o',
linewidth=2, markersize=23) # start point black edge
plt.plot(cell_list[i_sel][0][0],cell_list[i_sel][1][0], color = color_i, marker='o',
linewidth=2, markersize=21) # start point colored
plt.plot(cell_list[i_sel][0][0],cell_list[i_sel][1][0], color = 'w', marker='o',
linewidth=2, markersize=16) # start point white
plt.plot(cell_list[i_sel][0][0],cell_list[i_sel][1][0], color = 'k', marker=cell_marker_i,
linewidth=2, markersize=10*multiplying_factor) # start point number
#get current axes
ax = plt.gca()
#hide x-axis
ax.get_xaxis().set_visible(False)
#hide y-axis
ax.get_yaxis().set_visible(False)
# Save figure
#plt.savefig("Mosaic_TL_mCherry/Figure_in_vitro_tracking_mosaic_icons_lines_t1_more_white.png", dpi=600, bbox_inches='tight',pad_inches = 0)
# Plot figure
fig = plt.show()
# T = 31, 14 cells trajectories
%matplotlib inline
plt.rcParams['figure.figsize'] = [30, 20]
imgplot= mpimg.imread('Mosaic_TL_mCherry/MAX_prediction-t31-200epochsRGB-SB100um.png')
plt.imshow(imgplot)#, extent=[0, 300, 0, 300])
cmap = [#'#000000', #black
'#e69d00', #orange
'#f0e442', #yellow
'#009e73', #bluish green
'#56b4e9', #sky blue
'#cc79a7', #reddish purple
'#d55e00', #vermillion (reddish)
'#0072b2', #blue
#'#ffffff', #white
'#7f7f7f'] #grey
cell_marker = ['$1$','$2$','$3$','$4$','$5$','$6$','$7$','$8$','$9$','$10$','$11$','$12$','$13$','$14$']
n = 14
# extra lines
#border
x_points_lines = [cell_list[6][4][-1],cell_list[2][4][-1],cell_list[1][0][-1],cell_list[0][0][-1],
cell_list[5][0][-1],cell_list[11][0][-1],cell_list[13][0][-1],cell_list[12][0][-1],cell_list[9][0][-1],cell_list[6][4][-1]]
y_points_lines = [cell_list[6][5][-1],cell_list[2][5][-1],cell_list[1][1][-1],cell_list[0][1][-1],
cell_list[5][1][-1],cell_list[11][1][-1],cell_list[13][1][-1],cell_list[12][1][-1],cell_list[9][1][-1],cell_list[6][5][-1]]
#in between
x_points_lines_mid_1 = [cell_list[2][4][-1],cell_list[3][0][-1],cell_list[5][0][-1]]
y_points_lines_mid_1 = [cell_list[2][5][-1],cell_list[3][1][-1],cell_list[5][1][-1]]
x_points_lines_mid_2 = [cell_list[6][4][-1],cell_list[7][0][-1],cell_list[5][0][-1]]
y_points_lines_mid_2 = [cell_list[6][5][-1],cell_list[7][1][-1],cell_list[5][1][-1]]
x_points_lines_mid_3 = [cell_list[9][0][-1],cell_list[10][0][-1],cell_list[11][0][-1]]
y_points_lines_mid_3 = [cell_list[9][1][-1],cell_list[10][1][-1],cell_list[11][1][-1]]
plt.plot(x_points_lines, y_points_lines,'y',linewidth = 5)
plt.plot(x_points_lines_mid_1, y_points_lines_mid_1,'y',linewidth = 5)
plt.plot(x_points_lines_mid_2, y_points_lines_mid_2,'y',linewidth = 5)
plt.plot(x_points_lines_mid_3, y_points_lines_mid_3,'y',linewidth = 5)
for i in range(n):
if i >= 9:
multiplying_factor = 1.5 # use a larger icon for cells that start to count from 10
else:
multiplying_factor = 1
i_sel = i
cell_marker_i = cell_marker[i_sel]
if i < 8:
color_i = cmap[i]
if i >= 8 and i < 16:
color_i = cmap[i-8]
# end points
end_points = 1
if end_points == 1:
if cell_list[i_sel][2] == 'Dividing cell':
# daughter cell 1
plt.plot(cell_list[i_sel][4][-1],cell_list[i_sel][5][-1], color = 'k',
marker='D', linewidth=20, markersize=22) # end point black edge
plt.plot(cell_list[i_sel][4][-1],cell_list[i_sel][5][-1], color = color_i,
marker='D', linewidth=20, markersize=20) # end point colored
plt.plot(cell_list[i_sel][4][-1],cell_list[i_sel][5][-1], color = 'w',
marker='D', linewidth=20, markersize=15) # end point white
plt.plot(cell_list[i_sel][4][-1],cell_list[i_sel][5][-1], color = 'k', marker=cell_marker_i,
linewidth=20, markersize=8*multiplying_factor) # end point number
# daughter cell 2
plt.plot(cell_list[i_sel][6][-1],cell_list[i_sel][7][-1], color = 'k',
marker='D', linewidth=20, markersize=22) # end point black edge
plt.plot(cell_list[i_sel][6][-1],cell_list[i_sel][7][-1], color = color_i,
marker='D', linewidth=20, markersize=20) # end point colored
plt.plot(cell_list[i_sel][6][-1],cell_list[i_sel][7][-1], color = 'w',
marker='D', linewidth=20, markersize=15) # end point white
plt.plot(cell_list[i_sel][6][-1],cell_list[i_sel][7][-1], color = 'k', marker=cell_marker_i,
linewidth=20, markersize=10*multiplying_factor) # end point number
else: # no dividing cell
plt.plot(cell_list[i_sel][0][-1],cell_list[i_sel][1][-1], color = 'k',
marker='D', linewidth=20, markersize=22) # end point black edge
plt.plot(cell_list[i_sel][0][-1],cell_list[i_sel][1][-1], color = color_i,
marker='D', linewidth=20, markersize=20) # end point colored
plt.plot(cell_list[i_sel][0][-1],cell_list[i_sel][1][-1], color = 'w',
marker='D', linewidth=20, markersize=15) # end point white
plt.plot(cell_list[i_sel][0][-1],cell_list[i_sel][1][-1], color = 'k', marker=cell_marker_i,
linewidth=20, markersize=10*multiplying_factor) # end point number
#get current axes
ax = plt.gca()
#hide x-axis
ax.get_xaxis().set_visible(False)
#hide y-axis
ax.get_yaxis().set_visible(False)
# Save figure
#plt.savefig("Figure_in_vitro_tracking_mosaic_icons_lines_t31_endpointsonly.png", dpi=600, bbox_inches='tight',pad_inches = 0)
# Plot figure
fig = plt.show()
t_end_iv = 30 # t = 31
t_start_iv = 0 # t = 1
# Determine the length of the gastruloid at t_end_iv = 2 * long arm of ellipse fit. The long arm should be stored in the
# parameter 'b', and the short arm in parameter 'a', but the fitting tool sometimes stored them incorrectly.
# Therefore, check the value of a and b first to determine which parameter corresponds to the long arm.
# In our in vitro example, the long arm of the fitted ellipses should be on the y-axis. If a and b are incorrectly switched,
# the values stored in the x coordinates correspond to the long arm.
if a_iv[t_end_iv] > b_iv[t_end_iv]:
length_gastruloid_iv = 2*a_iv[t_end_iv] # used to normalize the distance the cells moved along the long axis
print('a > b, so we use the rotated x-coordinates to calculate the distance cells moved along the long axis of the gastruloid')
end_pos_on_long_axis_iv = rotated_cell_centers_iv[:,t_end_iv,0] # rotated x-coordinate at t = t_end_iv
end_pos_on_long_axis_iv_norm = (end_pos_on_long_axis_iv - xc_iv[t_end_iv])/(length_gastruloid_iv) # (rotated x-coordinate at t_end_iv - x_center_of_ellipse at t_end_iv) divided by length gastruloid
moved_end_iv_norm = ((end_pos_on_long_axis_iv - rotated_cell_centers_iv[:,t_start_iv,0]) - (xc_iv[t_end_iv] - xc_iv[t_start_iv]) ) / length_gastruloid_iv # to normalize, divide by length of gastruloid
else:
length_gastruloid_iv = 2*b_iv[t_end_iv] # used to normalize the distance the cells moved along the long axis
print('a < b, so we use the rotated y-coordinates to calculate the distance cells moved along the long axis of the gastruloid')
end_pos_on_long_axis_iv = rotated_cell_centers_iv[:,t_end_iv,1] # rotated y-coordinate at t = t_end_iv
end_pos_on_long_axis_iv_norm = (end_pos_on_long_axis_iv - yc_iv[t_end_iv])/(length_gastruloid_iv) # (rotated y-coordinate at t_end_iv - y_center_of_ellipse at t_end_iv) divided by length gastruloid
moved_end_iv_norm = ((end_pos_on_long_axis_iv - rotated_cell_centers_iv[:,t_start_iv,1]) - (yc_iv[t_end_iv] - yc_iv[t_start_iv]) ) / length_gastruloid_iv # to normalize, divide by length of gastruloid
# -------- FIGURE ----------
# Figure parameters
%matplotlib inline
font = {'size' : 12}
matplotlib.rc('font', **font)
matplotlib.rcParams['font.sans-serif'] = "Arial"
matplotlib.rcParams['font.family'] = "sans-serif"
plt.gcf().clear()
fig = plt.figure(1)
ax = fig.add_subplot(111)
plt.rcParams['figure.figsize'] = [10, 10]
plt.rcParams['xtick.labelsize'] = 'x-large'
plt.rcParams['ytick.labelsize'] = 'x-large'
plt.xlabel('End position projected on long axis', fontsize = 16)
plt.ylabel('Movement along long axis', fontsize = 16)
plt.xlim(-0.55, 0.55)
plt.ylim(-0.51, 0.51)
plt.xticks([-0.5, -0.25, 0, 0.25, 0.5],['-L/2','-L/4','0','L/4','L/2'])
plt.yticks([-0.5, -0.25, 0, 0.25, 0.5],['-L/2','-L/4','0','L/4','L/2'])
# In vitro time lapse data
ax.scatter(end_pos_on_long_axis_iv_norm, moved_end_iv_norm, color = '#cc79a7', marker = 'o', label = 'in vitro',s = 50)
handles, labels = ax.get_legend_handles_labels()
lgd = ax.legend(handles, labels, loc='upper left', bbox_to_anchor=(1,1))
fig.set_size_inches((110/5)/2.54, (76/5)/2.54)
plt.show()
a < b, so we use the rotated y-coordinates to calculate the distance cells moved along the long axis of the gastruloid
t_end_iv = 30 # t = 31
t_start_iv = 0 # t = 1
# Determine the length of the gastruloid at t_end_iv = 2 * long arm of ellipse fit. The long arm should be stored in the
# parameter 'b', and the short arm in parameter 'a', but the fitting tool sometimes stored them incorrectly.
# Therefore, check the value of a and b first to determine which parameter corresponds to the long or short arm.
# In our in vitro example, the long arm of the fitted ellipses should be on the y-axis. The short arm on the x-axis.
# If a and b are incorrectly switched, the values stored in the x coordinates correspond to the long arm, and y coordinates correspond to the short arm.
if a_iv[t_end_iv] > b_iv[t_end_iv]:
length_gastruloid_iv_short = 2*b_iv[t_end_iv] # used to normalize the distance the cells moved along the short axis
length_gastruloid_iv_long = 2*a_iv[t_end_iv] # used to normalize the position of the cells along the long axis
print('a > b, so we use the rotated y-coordinates to calculate the distance cells moved along the short axis of the gastruloid')
end_pos_on_short_axis_iv = rotated_cell_centers_iv[:,t_end_iv,1] # rotated y-coordinate at t = t_end_iv
end_pos_on_long_axis_iv = rotated_cell_centers_iv[:,t_end_iv,0] # rotated x-coordinate at t = t_end_iv
end_pos_on_short_axis_iv_norm = (end_pos_on_short_axis_iv - yc_iv[t_end_iv])/(length_gastruloid_iv_short) # (rotated y-coordinate at t_end_iv - y_center_of_ellipse at t_end_iv) divided by length gastruloid
end_pos_on_long_axis_iv_norm = (end_pos_on_long_axis_iv - xc_iv[t_end_iv])/(length_gastruloid_iv_long) # (rotated y-coordinate at t_end_iv - y_center_of_ellipse at t_end_iv) divided by length gastruloid
moved_end_iv_short_norm = ((end_pos_on_short_axis_iv - rotated_cell_centers_iv[:,t_start_iv,1]) - (yc_iv[t_end_iv] - yc_iv[t_start_iv]) ) / length_gastruloid_iv_short # to normalize, divide by length of gastruloid
moved_end_iv_long_norm = ((end_pos_on_long_axis_iv - rotated_cell_centers_iv[:,t_start_iv,0]) - (xc_iv[t_end_iv] - xc_iv[t_start_iv]) ) / length_gastruloid_iv_long # to normalize, divide by length of gastruloid
else:
length_gastruloid_iv_short = 2*a_iv[t_end_iv] # used to normalize the distance the cells moved along the short axis
length_gastruloid_iv_long = 2*b_iv[t_end_iv] # # used to normalize the position of the cells along the long axis
print('a < b, so we use the rotated x-coordinates to calculate the distance cells moved along the short axis of the gastruloid')
end_pos_on_short_axis_iv = rotated_cell_centers_iv[:,t_end_iv,0] # rotated x-coordinate at t = t_end_iv
end_pos_on_long_axis_iv = rotated_cell_centers_iv[:,t_end_iv,1] # rotated y-coordinate at t = t_end_iv
end_pos_on_short_axis_iv_norm = (end_pos_on_short_axis_iv - xc_iv[t_end_iv])/(length_gastruloid_iv_short) # (rotated x-coordinate at t_end_iv - x_center_of_ellipse at t_end_iv) divided by length gastruloid
end_pos_on_long_axis_iv_norm = (end_pos_on_long_axis_iv - yc_iv[t_end_iv])/(length_gastruloid_iv_long) # (rotated y-coordinate at t_end_iv - y_center_of_ellipse at t_end_iv) divided by length gastruloid
moved_end_iv_short_norm = ((end_pos_on_short_axis_iv - rotated_cell_centers_iv[:,t_start_iv,0]) - (xc_iv[t_end_iv] - xc_iv[t_start_iv]) ) / length_gastruloid_iv_short # to normalize, divide by length of gastruloid
moved_end_iv_long_norm = ((end_pos_on_long_axis_iv - rotated_cell_centers_iv[:,t_start_iv,1]) - (yc_iv[t_end_iv] - yc_iv[t_start_iv]) ) / length_gastruloid_iv_long # to normalize, divide by length of gastruloid
# -------- FIGURE ----------
# Figure parameters
%matplotlib inline
font = {'size' : 12}
matplotlib.rc('font', **font)
matplotlib.rcParams['font.sans-serif'] = "Arial"
matplotlib.rcParams['font.family'] = "sans-serif"
plt.gcf().clear()
fig = plt.figure(1)
ax = fig.add_subplot(111)
plt.rcParams['figure.figsize'] = [10, 10]
plt.rcParams['xtick.labelsize'] = 'x-large'
plt.rcParams['ytick.labelsize'] = 'x-large'
plt.xlabel('End position projected on short axis', fontsize = 16)
plt.ylabel('Movement along short axis', fontsize = 16)
plt.xlim(-1, 1)
plt.ylim(-1, 1)
plt.xticks([-1, -0.5, 0, 0.5, 1],['S','-S/2','0','S/2', 'S'])
plt.yticks([-1,-0.5, 0, 0.5, 1],['-S','-S/2','0','S/2','S'])
ax.scatter(end_pos_on_short_axis_iv_norm, moved_end_iv_short_norm, color = '#cc79a7', marker = 'o', s = 50)#label = 'in vitro',s = 50)
fig.set_size_inches((110/5)/2.54, (76/5)/2.54)
plt.show()
a < b, so we use the rotated x-coordinates to calculate the distance cells moved along the short axis of the gastruloid
# Read image
%matplotlib inline
img = cv2.imread('Mosaic_3/image0099900.png',0)
plt.imshow(img,cmap='gray')
plt.show()
# Erode the image multiple times, and subsequently dilate the image multiple times
%matplotlib inline
# Taking a matrix of size 5 as the kernel
kernel = np.ones((5, 5), np.uint8)
# The first parameter is the original image,
# kernel is the matrix with which image is
# convolved and third parameter is the number
# of iterations, which will determine how much
# you want to erode/dilate a given image.
img_erosion = cv2.erode(img, kernel, iterations=4)
img_dilation = cv2.dilate(img_erosion, kernel, iterations=4)
imagem = cv2.bitwise_not(img_erosion)
plt.imshow(imagem,cmap='gray')
plt.show()
plt.imshow(img_dilation,cmap='gray')
plt.show()
# save image
#cv2.imwrite('Mosaic_3/image0099900-binary.png', imagem)
# read the images and store them in a list in the correct order
image_path = Path('Mosaic_corrected')
#images = list(image_path.glob('*.png'))
image_list = []
for i in range(0,100000,100):
n = i
STRT = "Mosaic_corrected/image"
EXTN = ".png"
filename = STRT + str(n).zfill(7) + EXTN
image_list.append(imageio.imread(filename))
# from the list of images, write a GIF (standard 10 fps)
imageio.mimwrite('Mosaic_corrected/animated_mosaic_corrected.gif',image_list)
Alternatively, one coud also open the images in ImageJ, convert it to 8-bit, add a Monte Carlo step time stamp and save the sequence as a GIF.
def extract_outline():
for x in range(0,1000):
ind = list()
cont_all = list()
image_loc = 'Mosaic_corrected/image'+"{:05d}".format(x)+'00.png'
file_loc = 'Mosaic_corrected/outlines/outline_'+str(x) + '.txt'
#file1 = open(file_loc,"a")
mask = cv2.imread(image_loc)
mask = 255-mask
mask = mask[1:-1,1:-1]
mask = mask[1:-1,1
:-1]
mask = mask[1:-1,1:-1]
if mask is not None:
imgray = cv2.cvtColor(mask,cv2.COLOR_BGR2GRAY)
ret,thresh = cv2.threshold(imgray,200,255,0)
contours, im2 = cv2.findContours(thresh,cv2.RETR_TREE,cv2.CHAIN_APPROX_NONE)
maxsizeloc = 0
maxsize = 0
for i in range (len(contours)):
if contours[i].shape[0] > maxsize:
maxsizeloc = i
maxsize = contours[i].shape[0]
cont = np.zeros((contours[maxsizeloc].shape[0], 2))
cont[:,0] = contours[maxsizeloc][:,0,0]
cont[:,1] = contours[maxsizeloc][:,0,1]
plt.plot(cont[:,0], cont[:,1])
cont_all.append(cont)
#for i in range (maxsize):
#file1.write(str(cont_all[0][i][0]) + ',' + str(cont_all[0][i][1]) + '\n')
#file1.close()
return(cont_all)
extract_outline()
[array([[420., 296.],
[420., 297.],
[419., 298.],
...,
[423., 296.],
[422., 296.],
[421., 296.]])]
# Ellipse fit parameters
xc = np.empty(1000) # center coordinate x of ellipse
yc = np.empty(1000) # center coordinate y of ellipse
a = np.empty(1000) # length short arm of ellipse
b = np.empty(1000) # length long arm of ellipse
theta= np.empty(1000) # rotation angle of ellipse
# Figure parameters
%matplotlib inline
fig, axs = plt.subplots(figsize = (10,10))#, sharex=True, sharey=True)#, figsize=(5,15))
cmap = [#'#000000', #black
'#e69d00', #orange
'#56b4e9', #sky blue
'#009e73', #bluish green
'#f0e442', #yellow
'#0072b2', #blue
'#d55e00', #vermillion (reddish)
'#cc79a7', #reddish purple
#'#ffffff', #white
'#7f7f7f'] #grey
for i in range(1000):
# Load outline coordinates for MCS i
with open("Mosaic_corrected/outline_" + str(i) + ".txt") as file_name:
points = np.loadtxt(file_name, delimiter=",")
a_points = np.array(points)
x = a_points[:, 0]
y = a_points[:, 1]
#Calculate ellipse fit from outline and store parameters
ell = EllipseModel()
ell.estimate(a_points)
xc[i], yc[i], a[i], b[i], theta[i] = ell.params
# Write parameters to txt file
#file_loc = 'Mosaic_corrected/ellipse_params/MCS_'+str(i)+'-x_y_theta_a_b.txt'
#file1 = open(file_loc,"a")
#file1.write(str(xc[i]) + ',' + str(yc[i]) + ',' + str(theta[i]) + ',' + str(a[i]) + ',' + str(b[i]) + '\n')
# Plot ellipse fit for one MCS
t_end = 999
if i == t_end:
axs.scatter(x, y, color = cmap[4])
axs.set_ylim(0,1000)
axs.set_xlim(0,1000)
axs.scatter(xc[i], yc[i], color=cmap[0], s=50)
ell_patch = Ellipse((xc[i], yc[i]), 2*a[i], 2*b[i], theta[i]*180/np.pi, edgecolor=cmap[0], facecolor='none', linewidth = 3)
axs.add_patch(ell_patch)
if a[i] < b[i]:
print('a',a[i],' < b',b[i], 'theta = ', theta[i]*180/np.pi)
# Save figure
#plt.savefig("Mosaic_corrected/ellipse_fit_t999.png", dpi=600, pad_inches = 0)
plt.show()
print('Theta at MCS 100000 = ', theta[999]*180/np.pi)
Theta at MCS 100000 = 53.37818613761889
# Plot length of long arm and short arm over MCS - use to determine range of MCS to plot
%matplotlib inline
x1000 = list(range(0,1000))
for i in x1000:
if a[i] > b[i]: # a is long arm, b is short arm
plt.scatter(x1000[i],a[i], s = 0.5, color = 'k')
plt.scatter(x1000[i],b[i], s = 0.5, color = 'b')
elif a[i] < b[i]: # b is long arm, a is short arm
plt.scatter(x1000[i],b[i], s = 0.5, color = 'k')
plt.scatter(x1000[i],a[i], s = 0.5, color = 'b')
t_start = 199
line_1 = t_start
plt.plot([line_1, line_1],[0, 220])
plt.show()
# Define number of cells and the selected time points you want to plot the trajectory for
No_cells = 20 # number of cells that are being tracked
selected_points = [t_start,249,299,349,399,449,499,549,599,649,699,749,799,849,899,949,999] # 16 steps
# Define list of arrays for x- and y-coordinates for each tracked cell
path_x_sim = [None] * (No_cells)
path_y_sim = [None] * (No_cells)
path_x_sim_final = [None] * No_cells
path_y_sim_final = [None] * No_cells
# Store x- and y-coordinates for each tracked cell in lists 'path_x_sim_final' and 'path_y_sim_final'
for i in range(10,int(No_cells*10+10),10):
filepath = ""
file_dir = ['Mosaic_corrected/cell_','.txt']
filepath = file_dir[0]+str(i)+file_dir[1]
with open(filepath) as f:
for line in f:
currentline = line.split(",")
x_line = float(currentline[0])
y_line = float(currentline[1])
ii = int(i/10-1)
path_x_sim[ii] = np.append(path_x_sim[ii],x_line)
path_y_sim[ii] = np.append(path_y_sim[ii],y_line)
path_x_sim_final[ii] = path_x_sim[ii][1:(len(path_x_sim[ii])+1)]
path_y_sim_final[ii] = path_y_sim[ii][1:(len(path_x_sim[ii])+1)]
# Store x- and y-coordinates of each tracked cell for only the selected time points in lists 'x_sel_final' and 'y_sel_final'
x_sel = [None] * No_cells
y_sel = [None] * No_cells
x_sel_final = [None] * No_cells
y_sel_final = [None] * No_cells
for j in range(No_cells):
for k in range(len(selected_points)):
x_sel[j] = np.append(x_sel[j],2*path_x_sim_final[j][selected_points[k]])
y_sel[j] = np.append(y_sel[j],2*path_y_sim_final[j][selected_points[k]])
x_sel_final[j] = x_sel[j][1:(len(selected_points)+1)]
y_sel_final[j] = y_sel[j][1:(len(selected_points)+1)]
%matplotlib inline
# Select image (time point) and adjust size (optional)
imgplot_100= plt.imread('Mosaic_corrected/image0019900.png')
plt.rcParams['figure.figsize'] = [200, 100]
# Make image gray (so no dark blue cells)
R_100, G_100, B_100 = imgplot_100[:,:,0], imgplot_100[:,:,1], imgplot_100[:,:,2]
imgGray_100 = 0.1 * R_100 + 0.1 * G_100 + 0.7*B_100
cmap = [#'#000000', #black
'#e69d00', #orange
'#56b4e9', #sky blue
'#009e73', #bluish green
'#f0e442', #yellow
'#0072b2', #blue
'#d55e00', #vermillion (reddish)
'#7f7f7f', #grey
'#cc79a7']#, #reddish purple
cell_marker = ['$1$','$2$','$3$','$4$','$5$','$6$','$7$','$8$','$9$','$10$','$11$','$12$','$13$','$14$','$15$','$16$','$17$','$18$','$19$','$20$']
markers_numbered_selected = [11,10,14,17,16,18,3,5,6,7,8,9] # 12 cells
for i in range(len(markers_numbered_selected)):
if i <= 5:
multiplying_factor = 1.5 # use a larger icon for cells that start to count from 10
else:
multiplying_factor = 1
i_sel_rot = (markers_numbered_selected[i])*10-1
i_sel = markers_numbered_selected[i]-1
cell_marker_i = cell_marker[i_sel]
if i == 0:
color_i = cmap[2]
if i == 2:
color_i = cmap[7]
if i > 0 and i < 8 and i != 2:
color_i = cmap[i]
if i == 8:
color_i = cmap[0]
if i == 9:
color_i = cmap[2]
if i == 10:
color_i = cmap[1]
if i == 11:
color_i = cmap[3]
t_start_sel = 0
plt.plot(x_sel_final[i_sel][t_start_sel],y_sel_final[i_sel][t_start_sel], color = 'k', marker='o',
linewidth=2, markersize=130) # start point black edge
plt.plot(x_sel_final[i_sel][t_start_sel],y_sel_final[i_sel][t_start_sel], color = color_i, marker='o',
linewidth=2, markersize=120) # start point colored
plt.plot(x_sel_final[i_sel][t_start_sel],y_sel_final[i_sel][t_start_sel], color = 'w', marker='o',
linewidth=2, markersize=80) # start point white
plt.plot(x_sel_final[i_sel][t_start_sel],y_sel_final[i_sel][t_start_sel], color = 'k', marker=cell_marker_i,
linewidth=20, markersize=40*multiplying_factor) # start point number
# Save figure
plt.imshow(imgGray_100, cmap='gray')
# Save figure
#plt.savefig("Mosaic_corrected/Figure_in_silico_mosaic_corrected_12_icons_t199.png", bbox_inches='tight',pad_inches = 0)
plt.show()
%matplotlib inline
# Select image (time point) and adjust size (optional)
imgplot= plt.imread('Mosaic_corrected/image0099900.png') # final time point
plt.rcParams['figure.figsize'] = [200, 100]
# Make image gray (so no dark blue cells)
R, G, B = imgplot[:,:,0], imgplot[:,:,1], imgplot[:,:,2]
imgGray = 0.1 * R + 0.1 * G + 0.7*B #0.2989 * R + 0.5870 * G + 0.1140 * B
cmap = [#'#000000', #black
'#e69d00', #orange
'#56b4e9', #sky blue
'#009e73', #bluish green
'#f0e442', #yellow
'#0072b2', #blue
'#d55e00', #vermillion (reddish)
'#7f7f7f', #grey
'#cc79a7']#, #reddish purple
cell_marker = ['$1$','$2$','$3$','$4$','$5$','$6$','$7$','$8$','$9$','$10$','$11$','$12$','$13$','$14$','$15$','$16$','$17$','$18$','$19$','$20$']
markers_numbered_selected = [11,10,14,17,16,18,3,5,6,7,8,9]
for i in range(len(markers_numbered_selected)):
if i <= 5:
multiplying_factor = 1.5 # use a larger icon for cells that start to count from 10
else:
multiplying_factor = 1
i_sel_rot = (markers_numbered_selected[i])*10-1
i_sel = markers_numbered_selected[i]-1
cell_marker_i = cell_marker[i_sel]
if i == 0:
color_i = cmap[2]
if i == 2:
color_i = cmap[7]
if i > 0 and i < 8 and i != 2:
color_i = cmap[i]
if i == 8:
color_i = cmap[0]
if i == 9:
color_i = cmap[2]
if i == 10:
color_i = cmap[1]
if i == 11:
color_i = cmap[3]
plt.plot(x_sel_final[i_sel][-1],y_sel_final[i_sel][-1], color = 'k',
marker='D', linewidth=20, markersize=130) # end point black edge
plt.plot(x_sel_final[i_sel][-1],y_sel_final[i_sel][-1], color = color_i,
marker='D', linewidth=20, markersize=120) # end point colored
plt.plot(x_sel_final[i_sel][-1],y_sel_final[i_sel][-1], color = 'w',
marker='D', linewidth=20, markersize=70) # end point white
plt.plot(x_sel_final[i_sel][-1],y_sel_final[i_sel][-1], color = 'k', marker=cell_marker_i,
linewidth=20, markersize=40*multiplying_factor) # end point number
plt.imshow(imgGray, cmap='gray')
plt.show()
# Read in cell centers coordinates and store in cell_centers[cell][time][coordinate]
cell_centers = []
for i in range(1,201):
with open("Mosaic_corrected/cell_" + str(i) + ".txt") as file_name:
cell_centers.append(np.loadtxt(file_name, delimiter=","))
for cell in range(200):
for time in range(1000):
for coordinate in range(2):
cell_centers[cell][time][coordinate] = 2*cell_centers[cell][time][coordinate]
# Rotate cell centers and store in rotated_cell_centers (rotated around ellipse center)
rotated_cell_centers = np.empty((200, 1000, 4))
theta_corr = np.empty(len(theta))
for i in range(1,201): # cells
for j in range(len(cell_centers[0])): # time / MCS
# correct theta if it shows a jump (= orientation error)
if theta[j]*180/np.pi >64 : # gastruloid angle
theta_corr[j] = theta[j] - 1.57079633
elif theta[j]*180/np.pi <0 :
theta_corr[j] = theta[j] + 1.57079633
else:
theta_corr[j] = theta[j]
# Calculate rotated coordinates, rotate around center of ellips
# xc/yc = center coordinates of ellipse, O_x/O_y = coordinates of cell
O_x = np.cos(theta_corr[j]) * (cell_centers[i-1][j][0]-xc[j]) + np.sin(theta_corr[j]) * (cell_centers[i-1][j][1]-yc[j])
rotated_cell_centers[i-1,j,0] = O_x + xc[j] # rotated x-coordinate of cell
rotated_cell_centers[i-1,j,2] = O_x # x component of vector from center of ellipse to position of cell
O_y = np.cos(theta_corr[j]) * (cell_centers[i-1][j][1]-yc[j]) - np.sin(theta_corr[j]) * (cell_centers[i-1][j][0]-xc[j])
rotated_cell_centers[i-1,j,1] = O_y + yc[j] # rotated y-coordinate of cell
rotated_cell_centers[i-1,j,3] = O_y # y component of vector from center of ellipse to position of cell
# Investigate theta_corr
%matplotlib inline
print('Theta at t_start: ', theta[t_start]*180/np.pi, 'and theta at t_end: ',theta[999]*180/np.pi)
print('Corrected theta at t_start: ', theta_corr[t_start]*180/np.pi, 'and corrected theta at t_end: ',theta_corr[999]*180/np.pi)
theta_x = range(0,1000,1)
# Plot theta
plt.scatter(theta_x,theta*180/np.pi, s = 1, color = 'red')
# Plot theta_corr
plt.scatter(theta_x,theta_corr*180/np.pi, s = 1, color = 'black')
# Plot t_start
plt.scatter(199,theta_corr[199]*180/np.pi, s = 30, color = 'limegreen')
plt.xlim(-10, 1010)
Theta at t_start: 127.77361178960888 and theta at t_end: 53.37818613761889 Corrected theta at t_start: 37.773611605969975 and corrected theta at t_end: 53.37818613761889
(-10.0, 1010.0)
# CORRECTION FOR DRIFT AND FOR ROTATION
# Rotate rotated cell centers back to orientation gastruloid at final MCS, plot trajectories at final MCS and correct for drift
# Plot figure: each trajectory a different color (max. 8 colors), color in the color-blind-friendly palette
%matplotlib inline
# Select image (time point) and adjust size (optional)
imgplot= plt.imread('Mosaic_corrected/image0099900.png') # final time point
plt.rcParams['figure.figsize'] = [200, 100]
# Make image gray (so no dark blue cells)
R, G, B = imgplot[:,:,0], imgplot[:,:,1], imgplot[:,:,2]
imgGray = 0.1 * R + 0.1 * G + 0.7*B #0.2989 * R + 0.5870 * G + 0.1140 * B
# Calculate drift of gastruloid for each MCS w.r.t. final MCS using the difference in position of ellipse centers
x_corr_drift_all = np.empty((1000))
y_corr_drift_all = np.empty((1000))
# Rotate all rotated cell centers back to gastruloid's orientation at the final MCS (999)
rotated_cell_centers_rot = np.empty((200, 1000, 4))
cmap = [#'#000000', #black
'#e69d00', #orange
'#56b4e9', #sky blue
'#009e73', #bluish green
'#f0e442', #yellow
'#0072b2', #blue
'#d55e00', #vermillion (reddish)
'#7f7f7f', #grey
'#cc79a7']#, #reddish purple
#'#ffffff', #white
#'#7f7f7f'] #grey
cell_marker = ['$1$','$2$','$3$','$4$','$5$','$6$','$7$','$8$','$9$','$10$','$11$','$12$','$13$','$14$','$15$','$16$','$17$','$18$','$19$','$20$']
markers_numbered_selected = [11,10,14,17,16,18,3,5,6,7,8,9]
for i in range(0,200): # cells
for j in range(0,1000): # MCSs
x_corr_drift_all[j] = xc[-1]-xc[j]
y_corr_drift_all[j] = yc[-1]-yc[j]
O_x_rot = np.cos(2*np.pi - theta_corr[-1]) * (rotated_cell_centers[i][j][0]-xc[j]) + np.sin(2*np.pi - theta_corr[-1]) * (rotated_cell_centers[i][j][1]-yc[j])
rotated_cell_centers_rot[i,j,0] = O_x_rot + xc[j] # rotated x-coordinate of cell
rotated_cell_centers_rot[i,j,2] = O_x_rot # x component of vector from center of ellipse to position of cell
O_y_rot = np.cos(2*np.pi - theta_corr[-1]) * (rotated_cell_centers[i][j][1]-yc[j]) - np.sin(2*np.pi - theta_corr[-1]) * (rotated_cell_centers[i][j][0]-xc[j])
rotated_cell_centers_rot[i,j,1] = O_y_rot + yc[j] # rotated y-coordinate of cell
rotated_cell_centers_rot[i,j,3] = O_y_rot # y component of vector from center of ellipse to position of cell
for i in range(len(markers_numbered_selected)): # loop over cells
i_sel_rot = (markers_numbered_selected[i])*10-1
i_sel = (markers_numbered_selected[i])-1
cell_marker_i = cell_marker[i_sel]
if i == 0:
color_i = cmap[2]
if i == 2:
color_i = cmap[7]
if i > 0 and i < 8 and i != 2:
color_i = cmap[i]
if i == 8:
color_i = cmap[0]
if i == 9:
color_i = cmap[2]
if i == 10:
color_i = cmap[1]
if i == 11:
color_i = cmap[3]
plt.plot(rotated_cell_centers_rot[i_sel_rot,selected_points,0] + x_corr_drift_all[selected_points], rotated_cell_centers_rot[i_sel_rot,selected_points,1] + y_corr_drift_all[selected_points], color = 'k',
linewidth=25, markersize=25) # line connecting start and end point black
plt.plot(rotated_cell_centers_rot[i_sel_rot,selected_points,0] + x_corr_drift_all[selected_points], rotated_cell_centers_rot[i_sel_rot,selected_points,1] + y_corr_drift_all[selected_points], color = color_i,
linewidth=20, markersize=20) # line connecting start and end point, linewidth=25, markersize=25
for i in range(len(markers_numbered_selected)):
if i <= 5:
multiplying_factor = 1.5 # use a larger icon for cells that start to count from 10
else:
multiplying_factor = 1
i_sel_rot = (markers_numbered_selected[i])*10-1
i_sel = markers_numbered_selected[i]-1 #i
cell_marker_i = cell_marker[i_sel]
if i == 0:
color_i = cmap[2]
if i == 2:
color_i = cmap[7]
if i > 0 and i < 8 and i != 2:
color_i = cmap[i]
if i == 8:
color_i = cmap[0]
if i == 9:
color_i = cmap[2]
if i == 10:
color_i = cmap[1]
if i == 11:
color_i = cmap[3]
# Calculate drift for starting point w.r.t. the final point and add that to starting positions
x_corr_drift_f = xc[999]-xc[t_start]
y_corr_drift_f = yc[999]-yc[t_start]
# plot start point
plt.plot(rotated_cell_centers_rot[i_sel_rot,selected_points[0],0] + x_corr_drift_all[selected_points[0]], rotated_cell_centers_rot[i_sel_rot,selected_points[0],1] + y_corr_drift_all[selected_points[0]],
color = 'k', marker='o',linewidth=2, markersize=130) # start point black edge
plt.plot(rotated_cell_centers_rot[i_sel_rot,selected_points[0],0] + x_corr_drift_all[selected_points[0]], rotated_cell_centers_rot[i_sel_rot,selected_points[0],1] + y_corr_drift_all[selected_points[0]],
color = color_i, marker='o', linewidth=2, markersize=120) # start point colored, use alpha = .5 for transparency
plt.plot(rotated_cell_centers_rot[i_sel_rot,selected_points[0],0] + x_corr_drift_all[selected_points[0]], rotated_cell_centers_rot[i_sel_rot,selected_points[0],1] + y_corr_drift_all[selected_points[0]],
color = 'w', marker='o', linewidth=2, markersize=80)#70) # start point white
plt.plot(rotated_cell_centers_rot[i_sel_rot,selected_points[0],0] + x_corr_drift_all[selected_points[0]], rotated_cell_centers_rot[i_sel_rot,selected_points[0],1] + y_corr_drift_all[selected_points[0]],
color = 'k', marker=cell_marker_i, linewidth=20, markersize=40*multiplying_factor) # start point number
# plot end point
plt.plot(rotated_cell_centers_rot[i_sel_rot,selected_points[-1],0] + x_corr_drift_all[selected_points[-1]], rotated_cell_centers_rot[i_sel_rot,selected_points[-1],1] + y_corr_drift_all[selected_points[-1]],
color = 'k', marker='D', linewidth=20, markersize=130) # end point black edge
plt.plot(rotated_cell_centers_rot[i_sel_rot,selected_points[-1],0] + x_corr_drift_all[selected_points[-1]], rotated_cell_centers_rot[i_sel_rot,selected_points[-1],1] + y_corr_drift_all[selected_points[-1]],
color = color_i, marker='D', linewidth=20, markersize=120) # end point colored
plt.plot(rotated_cell_centers_rot[i_sel_rot,selected_points[-1],0] + x_corr_drift_all[selected_points[-1]], rotated_cell_centers_rot[i_sel_rot,selected_points[-1],1] + y_corr_drift_all[selected_points[-1]],
color = 'w', marker='D', linewidth=20, markersize=80)#70) # end point white
plt.plot(rotated_cell_centers_rot[i_sel_rot,selected_points[-1],0] + x_corr_drift_all[selected_points[-1]], rotated_cell_centers_rot[i_sel_rot,selected_points[-1],1] + y_corr_drift_all[selected_points[-1]],
color = 'k', marker=cell_marker_i,linewidth=20, markersize=40*multiplying_factor) # end point number
plt.imshow(imgGray, cmap='gray')
#plt.savefig("Mosaic_corrected/Figure_in_silico_mosaic_corrected_cdrift_crotation_12trajectories.png", bbox_inches='tight',pad_inches = 0)
plt.show()
# Define number of cells and the selected time points you want to plot the trajectory for
No_cells = 200 # number of cells that are being tracked
# Define list of arrays for x- and y-coordinates for each tracked cell
path_x_sim_all = [None] * (No_cells)
path_y_sim_all = [None] * (No_cells)
path_x_sim_all_final = [None] * No_cells
path_y_sim_all_final = [None] * No_cells
# Store x- and y-coordinates for each tracked cell in lists 'path_x_sim_all_final' and 'path_y_sim_all_final'
for i in range(1,int(No_cells)+1,1):
filepath = ""
file_dir = ['Mosaic_corrected/cell_','.txt']
filepath = file_dir[0]+str(i)+file_dir[1]
with open(filepath) as f:
for line in f:
currentline = line.split(",")
x_line = float(currentline[0])
y_line = float(currentline[1])
ii = int(i)-1
path_x_sim_all[ii] = np.append(path_x_sim_all[ii],x_line)
path_y_sim_all[ii] = np.append(path_y_sim_all[ii],y_line)
path_x_sim_all_final[ii] = path_x_sim_all[ii][1:(len(path_x_sim_all[ii])+1)]
path_y_sim_all_final[ii] = path_y_sim_all[ii][1:(len(path_y_sim_all[ii])+1)]
# Store x- and y-coordinates of each tracked cell for only the selected time points in lists 'x_sel_final' and 'y_sel_final'
x_sel_all = [None] * No_cells
y_sel_all = [None] * No_cells
x_sel_all_final = [None] * No_cells
y_sel_all_final = [None] * No_cells
for j in range(No_cells):
for k in range(len(selected_points)):
x_sel_all[j] = np.append(x_sel_all[j],2*path_x_sim_all_final[j][selected_points[k]]) # this coding can probably be improved
y_sel_all[j] = np.append(y_sel_all[j],2*path_y_sim_all_final[j][selected_points[k]])
x_sel_all_final[j] = x_sel_all[j][1:(len(selected_points)+1)]
y_sel_all_final[j] = y_sel_all[j][1:(len(selected_points)+1)]
# Color indicates orientation of direction vector
import matplotlib.colors as colors
import matplotlib.cm as cmx
import matplotlib as mpl
name = 'twilight'# or 'hsv'
cmap=mpl.colormaps[name]
#cmap = plt.cm.jet
cNorm = colors.Normalize(vmin=-np.pi, vmax=np.pi)
scalarMap = cmx.ScalarMappable(norm=cNorm,cmap=cmap)
%matplotlib inline
imgplot_999= plt.imread('Mosaic_corrected/image0099900.png')
plt.rcParams['figure.figsize'] = [200, 100]
# Make image gray (so no dark blue cells)
R_999, G_999, B_999 = imgplot_999[:,:,0], imgplot_999[:,:,1], imgplot_999[:,:,2]
imgGray_999 = 0.1 * R_999 + 0.1 * G_999 + 0.7*B_999
for i in range(No_cells):
# Rotated and rotated back, vector downscaled by factor 20
x_dx_rot_rot_drift = ((rotated_cell_centers_rot[i,t_end,0] - rotated_cell_centers_rot[i,t_start,0]) - (xc[t_end] - xc[t_start]))/20
y_dy_rot_rot_drift = ((rotated_cell_centers_rot[i,t_end,1] - rotated_cell_centers_rot[i,t_start,1]) - (yc[t_end] - yc[t_start]))/20
x_arrow_base_rot_rot = rotated_cell_centers_rot[i,t_end,0] - x_dx_rot_rot_drift
y_arrow_base_rot_rot = rotated_cell_centers_rot[i,t_end,1] - y_dy_rot_rot_drift
angle_cell = np.arctan2((y_dy_rot_rot_drift),(x_dx_rot_rot_drift))
colorVal = scalarMap.to_rgba(angle_cell)
plt.arrow(x_arrow_base_rot_rot, y_arrow_base_rot_rot, x_dx_rot_rot_drift, y_dy_rot_rot_drift, width = 1, head_width=6, head_length=6, color=colorVal)
plt.arrow(x_arrow_base_rot_rot+300, y_arrow_base_rot_rot-290, x_dx_rot_rot_drift, y_dy_rot_rot_drift, width = 1, head_width=6, head_length=6, color=colorVal)
# Save figure
plt.imshow(imgGray_999, cmap='gray', alpha = .5)
#plt.savefig("Mosaic_corrected/Figure_in_silico_mosaic_corrected_t999_all_vectors_corr_drift_rot_DIV20_angle_twilight.png", bbox_inches='tight',pad_inches = 0)
plt.show()
# https://stackoverflow.com/questions/18748328/plotting-arrows-with-different-color-in-matplotlib
# use https://stackoverflow.com/questions/19576495/color-matplotlib-quiver-field-according-to-magnitude-and-direction
# Plot color wheel
%matplotlib inline
from matplotlib import cm
import matplotlib as mpl
fig = plt.figure()
display_axes = fig.add_axes([0.1,0.1,0.8,0.8], projection='polar')
norm = mpl.colors.Normalize(0,2*np.pi)
quant_steps = 2056
cb = mpl.colorbar.ColorbarBase(display_axes, cmap=cm.get_cmap('twilight_shifted',quant_steps),norm=norm, orientation='horizontal')
# aesthetics - get rid of border and axis labels
cb.outline.set_visible(False)
display_axes.set_axis_off()
plt.savefig("Mosaic_corrected/color_wheel_twilight.png", bbox_inches='tight',pad_inches = 0)
plt.show()
/var/folders/9w/72zthbm57f33y8wq5zqg5kjw0000gn/T/ipykernel_1167/1637029653.py:12: MatplotlibDeprecationWarning: Auto-removal of grids by pcolor() and pcolormesh() is deprecated since 3.5 and will be removed two minor releases later; please call grid(False) first.
cb = mpl.colorbar.ColorbarBase(display_axes, cmap=cm.get_cmap('twilight_shifted',quant_steps),norm=norm, orientation='horizontal')
# Determine the length of the simulated gastruloid at t_end = 2 * long arm of ellipse fit. The long arm should be stored in the
# parameter 'b', and the short arm in parameter 'a', but the fitting tool sometimes stored them incorrectly.
# Therefore, check the value of a and b first to determine which parameter corresponds to the long or short arm.
# In our in silico example, the long arm of the fitted ellipses is on the x-axis. The short arm on the y-axis.
# If a and b are incorrectly switched, the values stored in the y coordinates correspond to the long arm, and x coordinates correspond to the short arm.
if a[t_end] > b[t_end]:
print('a > b, so x-axis is long axis')
length_gastruloid = 2*a[t_end] # used to normalize the position of the cells along the long axis
end_pos_on_long_axis = rotated_cell_centers[:,t_end,0] # rotated x-coordinate at t = t_end
end_pos_on_long_axis_norm = (end_pos_on_long_axis - xc[t_end])/(length_gastruloid) # (rotated x-coordinate at t_end - x_center_of_ellipse at t_end) divided by length gastruloid
# (rotated x-coordinate at t_end - rotated x-coordinate at t_start) - difference between center of ellipses to correct for drift
moved_end_norm = ((end_pos_on_long_axis - rotated_cell_centers[:,t_start,0]) - (xc[t_end] - xc[t_start]) ) / length_gastruloid # divided by length of gastruloid to normalize
else:
print('a < b, so y-axis is long axis')
length_gastruloid = 2*b[t_end]
# Rotated cell centers around ellipse center
end_pos_on_long_axis = rotated_cell_centers[:,t_end,1] # rotated y-coordinate at t = t_end
end_pos_on_long_axis_norm = (end_pos_on_long_axis - yc[t_end])/(length_gastruloid) # (rotated y-coordinate at t_end - y_center_of_ellipse at t_end) divided by length gastruloid
# rotated y-coordinate at t_end - rotated y-coordinate at t_start
moved_end_norm = ((end_pos_on_long_axis - rotated_cell_centers[:,t_start,1]) - (yc[t_end] - yc[t_start]) ) / length_gastruloid # to normalize, divide by length of gastruloid
# In vitro time lapse data
# See cells under "In vitro"
# -------- FIGURE ----------
# Figure parameters
%matplotlib inline
font = {'size' : 12}
matplotlib.rc('font', **font)
matplotlib.rcParams['font.sans-serif'] = "Arial"
matplotlib.rcParams['font.family'] = "sans-serif"
plt.gcf().clear()
fig = plt.figure(1)
ax = fig.add_subplot(111)
plt.rcParams['figure.figsize'] = [10, 10]
plt.rcParams['xtick.labelsize'] = 'x-large'
plt.rcParams['ytick.labelsize'] = 'x-large'
plt.xlabel('End position projected on long axis', fontsize = 16)
plt.ylabel('Movement along long axis', fontsize = 16)
plt.xlim(-0.55, 0.55)
plt.ylim(-0.51, 0.51)
plt.xticks([-0.5, -0.25, 0, 0.25, 0.5],['-L/2','-L/4','0','L/4','L/2'])
plt.yticks([-0.5, -0.25, 0, 0.25, 0.5],['-L/2','-L/4','0','L/4','L/2'])
# In vitro time lapse data
ax.scatter(end_pos_on_long_axis_iv_norm, moved_end_iv_norm, color = '#cc79a7', marker = 'o', label = 'in vitro',s = 50)
# Plot rotated cell centers around ellipse center
ax.scatter(end_pos_on_long_axis_norm, moved_end_norm, s = 5, color = 'k', label = 'in silico')
handles, labels = ax.get_legend_handles_labels()
lgd = ax.legend(handles, labels, loc='upper left', bbox_to_anchor=(1,1))
fig.set_size_inches((110/5)/2.54, (76/5)/2.54)
# Save figure
#plt.savefig("Mosaic_corrected/Figure_corr_movement_along_long_axis_legend.svg", bbox_extra_artists=([lgd]), bbox_inches='tight')
plt.show()
a > b, so x-axis is long axis
# Plot movement along short axis
if a[t_end] > b[t_end]: # if a is long axis, rotated y-axis is short axis
print('a > b, so y-axis is short axis')
length_gastruloid_short = 2*b[t_end] # short axis
# Rotated cell centers around ellipse center
end_pos_on_short_axis = rotated_cell_centers[:,t_end,1] # rotated y-coordinate at t = t_end
end_pos_on_short_axis_norm = (end_pos_on_short_axis - yc[t_end])/(length_gastruloid_short) # (rotated y-coordinate at t_end - y_center_of_ellipse at t_end) divided by (short) length gastruloid
# (rotated y-coordinate at t_end - rotated y-coordinate at t_start) - difference between center of ellipses to correct for drift
moved_end_short_norm = ((end_pos_on_short_axis - rotated_cell_centers[:,t_start,1]) - (yc[t_end] - yc[t_start]) ) / length_gastruloid_short # divided by length of gastruloid to normalize
else: # if b is long axis, rotated x-axis is short axis
print('a < b, so x-axis is short axis')
length_gastruloid_short = 2*a[t_end]
# Rotated cell centers around ellipse center
end_pos_on_short_axis = rotated_cell_centers[:,t_end,0] # rotated x-coordinate at t = t_end
end_pos_on_short_axis_norm = (end_pos_on_short_axis - xc[t_end])/(length_gastruloid_short) # (rotated y-coordinate at t_end - y_center_of_ellipse at t_end) divided by length gastruloid
# rotated x-coordinate at t_end - rotated x-coordinate at t_start
moved_end_short_norm = ((end_pos_on_short_axis - rotated_cell_centers[:,t_start,0]) - (xc[t_end] - xc[t_start]) ) / length_gastruloid_short # to normalize, divide by length of gastruloid
# In vitro time lapse data
# See cells under "In vitro"
# -------- FIGURE ----------
# Figure parameters
%matplotlib inline
font = {'size' : 12}
matplotlib.rc('font', **font)
matplotlib.rcParams['font.sans-serif'] = "Arial"
matplotlib.rcParams['font.family'] = "sans-serif"
plt.gcf().clear()
fig = plt.figure(1)
ax = fig.add_subplot(111)
plt.rcParams['figure.figsize'] = [10, 10]
plt.rcParams['xtick.labelsize'] = 'x-large'
plt.rcParams['ytick.labelsize'] = 'x-large'
plt.xlabel('End position projected on short axis', fontsize = 16)
plt.ylabel('Movement along short axis', fontsize = 16)
plt.xlim(-1, 1)
plt.ylim(-1, 1)
plt.xticks([-1, -0.5, 0, 0.5, 1],['S','-S/2','0','S/2', 'S'])
plt.yticks([-1,-0.5, 0, 0.5, 1],['-S','-S/2','0','S/2','S'])
ax.scatter(end_pos_on_short_axis_iv_norm, moved_end_iv_short_norm, color = '#cc79a7', marker = 'o', s = 50)#label = 'in vitro',s = 50)
ax.scatter(end_pos_on_short_axis_norm, moved_end_short_norm, s = 5, color = 'k')#, label = 'in silico')
fig.set_size_inches((110/5)/2.54, (76/5)/2.54)
# Save figure
#plt.savefig("Mosaic_corrected/Figure_corr_movement_along_short_axis_end_pos_short_eqaxS.svg", bbox_extra_artists=([lgd]), bbox_inches='tight')
plt.show()
a > b, so y-axis is short axis
# Plot movement along short axis
if a[t_end] > b[t_end]: # if a is long axis, rotated y-axis is short axis
print('a > b, so y-axis is short axis')
length_gastruloid_short = 2*b[t_end] # short axis
# Rotated cell centers around ellipse center
end_pos_on_short_axis = rotated_cell_centers[:,t_end,1] # rotated y-coordinate at t = t_end
end_pos_on_short_axis_norm = (end_pos_on_short_axis - yc[t_end])/(length_gastruloid_short) # (rotated y-coordinate at t_end - y_center_of_ellipse at t_end) divided by (short) length gastruloid
# (rotated y-coordinate at t_end - rotated y-coordinate at t_start) - difference between center of ellipses to correct for drift
moved_end_short_norm = ((end_pos_on_short_axis - rotated_cell_centers[:,t_start,1]) - (yc[t_end] - yc[t_start]) ) / length_gastruloid_short # divided by length of gastruloid to normalize
else: # if b is long axis, rotated x-axis is short axis
print('a < b, so x-axis is short axis')
length_gastruloid_short = 2*a[t_end]
# Rotated cell centers around ellipse center
end_pos_on_short_axis = rotated_cell_centers[:,t_end,0] # rotated x-coordinate at t = t_end
end_pos_on_short_axis_norm = (end_pos_on_short_axis - xc[t_end])/(length_gastruloid_short) # (rotated y-coordinate at t_end - y_center_of_ellipse at t_end) divided by length gastruloid
# rotated x-coordinate at t_end - rotated x-coordinate at t_start
moved_end_short_norm = ((end_pos_on_short_axis - rotated_cell_centers[:,t_start,0]) - (xc[t_end] - xc[t_start]) ) / length_gastruloid_short # to normalize, divide by length of gastruloid
inw_outw = np.arange(0)
for i in range(0,len(end_pos_on_short_axis_norm)):
#print(i)
if end_pos_on_short_axis_norm[i] > 0 and moved_end_short_norm[i] < 0:
inw_outw = np.append(inw_outw, 1) # inwards
elif end_pos_on_short_axis_norm[i] > 0 and moved_end_short_norm[i] > 0 and moved_end_short_norm[i] > end_pos_on_short_axis_norm[i]:
inw_outw = np.append(inw_outw, 2) # inwards, outwards
elif end_pos_on_short_axis_norm[i] > 0 and moved_end_short_norm[i] > 0 and moved_end_short_norm[i] < end_pos_on_short_axis_norm[i]:
inw_outw = np.append(inw_outw, 0) # outwards
elif end_pos_on_short_axis_norm[i] < 0 and moved_end_short_norm[i] < 0 and moved_end_short_norm[i] < end_pos_on_short_axis_norm[i]:
inw_outw = np.append(inw_outw, 2) # inwards, outwards
elif end_pos_on_short_axis_norm[i] < 0 and moved_end_short_norm[i] < 0 and moved_end_short_norm[i] > end_pos_on_short_axis_norm[i]:
inw_outw = np.append(inw_outw, 0) # outwards
elif end_pos_on_short_axis_norm[i] < 0 and moved_end_short_norm[i] > 0:
inw_outw = np.append(inw_outw, 1) # inwards
else:
inw_outw = np.append(inw_outw, 3) # something is zero
# In vitro time lapse data
# See cells under "In vitro"
# -------- FIGURE ----------
# Figure parameters
%matplotlib inline
font = {'size' : 12}
matplotlib.rc('font', **font)
matplotlib.rcParams['font.sans-serif'] = "Arial"
matplotlib.rcParams['font.family'] = "sans-serif"
plt.gcf().clear()
fig = plt.figure(1)
ax = fig.add_subplot(111)
plt.rcParams['figure.figsize'] = [10, 10]
plt.rcParams['xtick.labelsize'] = 'x-large'
plt.rcParams['ytick.labelsize'] = 'x-large'
plt.xlabel('End position projected on short axis', fontsize = 16)
plt.ylabel('Movement along short axis', fontsize = 16)
plt.xlim(-1, 1)
plt.ylim(-1, 1)
plt.xticks([-1, -0.5, 0, 0.5, 1],['S','-S/2','0','S/2', 'S'])
plt.yticks([-1,-0.5, 0, 0.5, 1],['-S','-S/2','0','S/2','S'])
# In vitro time lapse data
inw_outw_iv = np.arange(0)
for i in range(0,len(end_pos_on_short_axis_iv_norm)):
#print(i)
if end_pos_on_short_axis_iv_norm[i] > 0 and moved_end_iv_short_norm[i] < 0:
inw_outw_iv = np.append(inw_outw_iv, 1) # inwards
elif end_pos_on_short_axis_iv_norm[i] > 0 and moved_end_iv_short_norm[i] > 0 and moved_end_iv_short_norm[i] > end_pos_on_short_axis_iv_norm[i]:
inw_outw_iv = np.append(inw_outw_iv, 2) # inwards, outwards
elif end_pos_on_short_axis_iv_norm[i] > 0 and moved_end_iv_short_norm[i] > 0 and moved_end_iv_short_norm[i] < end_pos_on_short_axis_iv_norm[i]:
inw_outw_iv = np.append(inw_outw_iv, 0) # outwards
elif end_pos_on_short_axis_iv_norm[i] < 0 and moved_end_iv_short_norm[i] < 0 and moved_end_iv_short_norm[i] < end_pos_on_short_axis_iv_norm[i]:
inw_outw_iv = np.append(inw_outw_iv, 2) # inwards, outwards
elif end_pos_on_short_axis_iv_norm[i] < 0 and moved_end_iv_short_norm[i] < 0 and moved_end_iv_short_norm[i] > end_pos_on_short_axis_iv_norm[i]:
inw_outw_iv = np.append(inw_outw_iv, 0) # outwards
elif end_pos_on_short_axis_iv_norm[i] < 0 and moved_end_iv_short_norm[i] > 0:
inw_outw_iv = np.append(inw_outw_iv, 1) # inwards
else:
inw_outw_iv = np.append(inw_outw_iv, 3) # something is zero
color_i = ['blue', #outwards
'orange', #inwards
'skyblue', #inwards outwards
'yellow'] # not sure]
for i in range(0,len(end_pos_on_long_axis_iv_norm)):
ax.scatter(end_pos_on_short_axis_iv_norm[i], moved_end_iv_short_norm[i], color = color_i[inw_outw_iv[i]], marker = 'o',s = 50)
for i in range(0,len(end_pos_on_short_axis_norm)):
ax.scatter(end_pos_on_short_axis_norm[i], moved_end_short_norm[i], s = 5, color = color_i[inw_outw[i]])#, label = 'in silico')
fig.set_size_inches((110/5)/2.54, (76/5)/2.54)
# Save figure
#plt.savefig("Mosaic_corrected/Figure_corr_movement_along_short_axis_end_pos_short_colored_eqaxS.svg", bbox_extra_artists=([lgd]), bbox_inches='tight')
plt.show()
a > b, so y-axis is short axis
# Plot movement along short axis vs. end position projected on long axis
if a[t_end] > b[t_end]: # if a is long axis, rotated y-axis is short axis
print('a > b, so y-axis is short axis')
length_gastruloid_short = 2*b[t_end] # short axis
length_gastruloid_long = 2*a[t_end]
# Rotated cell centers around ellipse center
end_pos_on_short_axis = rotated_cell_centers[:,t_end,1] # rotated y-coordinate at t = t_end
end_pos_on_long_axis = rotated_cell_centers[:,t_end,0] # rotated x-coordinate at t = t_end
end_pos_on_short_axis_norm = (end_pos_on_short_axis - yc[t_end])/(length_gastruloid_short) # (rotated y-coordinate at t_end - y_center_of_ellipse at t_end) divided by (short) length gastruloid
end_pos_on_long_axis_norm = (end_pos_on_long_axis - xc[t_end])/(length_gastruloid_long) # (rotated x-coordinate at t_end - y_center_of_ellipse at t_end) divided by (short) length gastruloid
# (rotated y-coordinate at t_end - rotated y-coordinate at t_start) - difference between center of ellipses to correct for drift
moved_end_short_norm = ((end_pos_on_short_axis - rotated_cell_centers[:,t_start,1]) - (yc[t_end] - yc[t_start]) ) / length_gastruloid_short # divided by length of gastruloid to normalize
moved_end_long_norm = ((end_pos_on_long_axis - rotated_cell_centers[:,t_start,0]) - (xc[t_end] - xc[t_start]) ) / length_gastruloid_long # divided by length of gastruloid to normalize
else: # if b is long axis, rotated x-axis is short axis
print('a < b, so x-axis is short axis, y-axis is long axis')
length_gastruloid_short = 2*a[t_end]
length_gastruloid_long = 2*b[t_end]
# Rotated cell centers around ellipse center
end_pos_on_short_axis = rotated_cell_centers[:,t_end,0] # rotated x-coordinate at t = t_end
end_pos_on_long_axis = rotated_cell_centers[:,t_end,1] # rotated y-coordinate at t = t_end
end_pos_on_short_axis_norm = (end_pos_on_short_axis - xc[t_end])/(length_gastruloid_short) # (rotated y-coordinate at t_end - y_center_of_ellipse at t_end) divided by length gastruloid
end_pos_on_long_axis_norm = (end_pos_on_long_axis - yc[t_end])/(length_gastruloid_long) # (rotated x-coordinate at t_end - y_center_of_ellipse at t_end) divided by (short) length gastruloid
# rotated x-coordinate at t_end - rotated x-coordinate at t_start
moved_end_short_norm = ((end_pos_on_short_axis - rotated_cell_centers[:,t_start,0]) - (xc[t_end] - xc[t_start]) ) / length_gastruloid_short # to normalize, divide by length of gastruloid
moved_end_long_norm = ((end_pos_on_long_axis - rotated_cell_centers[:,t_start,1]) - (yc[t_end] - yc[t_start]) ) / length_gastruloid_long # divided by length of gastruloid to normalize
inw_outw = np.arange(0)
for i in range(0,len(end_pos_on_short_axis_norm)):
#print(i)
if end_pos_on_short_axis_norm[i] > 0 and moved_end_short_norm[i] < 0:
inw_outw = np.append(inw_outw, 1) # inwards
elif end_pos_on_short_axis_norm[i] > 0 and moved_end_short_norm[i] > 0 and moved_end_short_norm[i] > end_pos_on_short_axis_norm[i]:
inw_outw = np.append(inw_outw, 2) # inwards, outwards
elif end_pos_on_short_axis_norm[i] > 0 and moved_end_short_norm[i] > 0 and moved_end_short_norm[i] < end_pos_on_short_axis_norm[i]:
inw_outw = np.append(inw_outw, 0) # outwards
elif end_pos_on_short_axis_norm[i] < 0 and moved_end_short_norm[i] < 0 and moved_end_short_norm[i] < end_pos_on_short_axis_norm[i]:
inw_outw = np.append(inw_outw, 2) # inwards, outwards
elif end_pos_on_short_axis_norm[i] < 0 and moved_end_short_norm[i] < 0 and moved_end_short_norm[i] > end_pos_on_short_axis_norm[i]:
inw_outw = np.append(inw_outw, 0) # outwards
elif end_pos_on_short_axis_norm[i] < 0 and moved_end_short_norm[i] > 0:
inw_outw = np.append(inw_outw, 1) # inwards
else:
inw_outw = np.append(inw_outw, 3) # something is zero
# In vitro time lapse data
# See cells under "In vitro"
inw_outw_iv = np.arange(0)
for i in range(0,len(end_pos_on_short_axis_iv_norm)):
#print(i)
if end_pos_on_short_axis_iv_norm[i] > 0 and moved_end_iv_short_norm[i] < 0:
inw_outw_iv = np.append(inw_outw_iv, 1) # inwards
elif end_pos_on_short_axis_iv_norm[i] > 0 and moved_end_iv_short_norm[i] > 0 and moved_end_iv_short_norm[i] > end_pos_on_short_axis_iv_norm[i]:
inw_outw_iv = np.append(inw_outw_iv, 2) # inwards, outwards
elif end_pos_on_short_axis_iv_norm[i] > 0 and moved_end_iv_short_norm[i] > 0 and moved_end_iv_short_norm[i] < end_pos_on_short_axis_iv_norm[i]:
inw_outw_iv = np.append(inw_outw_iv, 0) # outwards
elif end_pos_on_short_axis_iv_norm[i] < 0 and moved_end_iv_short_norm[i] < 0 and moved_end_iv_short_norm[i] < end_pos_on_short_axis_iv_norm[i]:
inw_outw_iv = np.append(inw_outw_iv, 2) # inwards, outwards
elif end_pos_on_short_axis_iv_norm[i] < 0 and moved_end_iv_short_norm[i] < 0 and moved_end_iv_short_norm[i] > end_pos_on_short_axis_iv_norm[i]:
inw_outw_iv = np.append(inw_outw_iv, 0) # outwards
elif end_pos_on_short_axis_iv_norm[i] < 0 and moved_end_iv_short_norm[i] > 0:
inw_outw_iv = np.append(inw_outw_iv, 1) # inwards
else:
inw_outw_iv = np.append(inw_outw_iv, 3) # something is zero
# -------- FIGURE ----------
# Figure parameters
%matplotlib inline
font = {'size' : 12}
matplotlib.rc('font', **font)
matplotlib.rcParams['font.sans-serif'] = "Arial"
matplotlib.rcParams['font.family'] = "sans-serif"
plt.gcf().clear()
fig = plt.figure(1)
ax = fig.add_subplot(111)
plt.rcParams['figure.figsize'] = [10, 10]
plt.rcParams['xtick.labelsize'] = 'x-large'
plt.rcParams['ytick.labelsize'] = 'x-large'
plt.xlabel('End position projected on long axis', fontsize = 16)
plt.ylabel('Movement along short axis', fontsize = 16)
plt.xlim(-0.55, 0.55)
plt.ylim(-1, 1)#-0.38, 0.38)
plt.xticks([-0.5, -0.25, 0, 0.25, 0.5],['-L/2','-L/4','0','L/4','L/2'])
plt.yticks([-1,-0.5, 0, 0.5,1],['-S','-S/2','0','S/2','S'])
color_i = ['blue', #outwards
'orange', #inwards
'skyblue', #inwards outwards
'yellow'] # not sure]
# In vitro time lapse data
for i in range(0,len(end_pos_on_long_axis_iv_norm)):
ax.scatter(end_pos_on_long_axis_iv_norm[i], moved_end_iv_short_norm[i], color = color_i[inw_outw_iv[i]], marker = 'o',s = 50)
# Plot rotated cell centers around ellipse center
for i in range(0,len(end_pos_on_short_axis_norm)):
ax.scatter(end_pos_on_long_axis_norm[i], moved_end_short_norm[i], s = 5, color = color_i[inw_outw[i]])#, label = 'in silico')
fig.set_size_inches((110/5)/2.54, (76/5)/2.54)
# Save figure
#plt.savefig("Mosaic_corrected/Figure_corr_movement_along_short_axis_end_pos_long_colored_more_spaceS.svg", bbox_extra_artists=([lgd]), bbox_inches='tight')
plt.show()
a > b, so y-axis is short axis
# First, change image type from 16 bit (raw tiff file) or 32 bit (prediction) to 8 bit, e.g. in FIJI.
# Read RGB images
%matplotlib inline
# Read raw and predicted image
#image_raw = cv2.imread('Mosaic_TL_mCherry/N2V_comp/T1-z11-raw-crop_8bit.tif',0)
#image_pred = cv2.imread('Mosaic_TL_mCherry/N2V_comp/prediction-t1-200epochs-z11_8bit_crop.tif',0)
#image_raw = cv2.imread('Mosaic_TL_mCherry/N2V_comp/t5-z10-crop_8bit.tif',0)
#image_pred = cv2.imread('Mosaic_TL_mCherry/N2V_comp/prediction-t5_t7-model-200epochs-z10-crop_8bit.tif',0)
#image_raw = cv2.imread('Mosaic_TL_mCherry/N2V_comp/t23-z6-crop_8bit.tif',0)
#image_pred = cv2.imread('Mosaic_TL_mCherry/N2V_comp/prediction-t23_t31-model-200epochs_z6_crop_8bit.tif',0)
image_raw = cv2.imread('Mosaic_TL_mCherry/N2V_comp/MAX_T31_crop-raw_8bit.tif',0)
image_pred = cv2.imread('Mosaic_TL_mCherry/N2V_comp/MAX_prediction-t31-200epochs-crop_8bit.tif',0)
plt.imshow(image_raw,cmap='gray')
plt.show()
plt.imshow(image_pred,cmap='gray')
plt.show()
# create the histogram
counts_raw, bins_raw = np.histogram(image_raw, bins=256, range=(0, 255))
counts_pred, bins_pred = np.histogram(image_pred, bins=256, range=(0, 255))
%matplotlib inline
plt.bar(bins_raw[:-1],counts_raw,width=1)
plt.bar(bins_pred[:-1],counts_pred,width=1)
<BarContainer object of 256 artists>
# flatten the 2D image arrays
image_raw_flat = image_raw.flatten()
image_pred_flat = image_pred.flatten()
# Calculate 5th percentile of raw data, such that that value can be used as the minimum displayed value
image_raw_sort = np.sort(image_raw_flat)
print(image_raw_sort)
q5_raw=int(np.round(0.05*(len(image_raw_flat))))
print(q5_raw)
print('minimum displayed value for raw 8-bit image is: ', image_raw_sort[q5_raw])
# Calculate 5th percentile of prediction data, such that that value can be used as the minimum displayed value
image_pred_sort = np.sort(image_pred_flat)
print(image_pred_sort)
q5_pred=int(np.round(0.05*(len(image_pred_flat))))
print(q5_pred)
print('minimum displayed value for pred. 8-bit image is: ', image_pred_sort[q5_pred])
[ 2 3 3 ... 255 255 255] 22714 minimum displayed value for raw 8-bit image is: 16 [ 1 1 1 ... 255 255 255] 22714 minimum displayed value for pred. 8-bit image is: 4